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1. 


INTRODUCTION 


This  is  the  final  report  for  the  project  entitled  "Computer 
Tomography  and  Hybrid  Optical/Digital  Methods  for  Aerodynamic 
Measurements"  sponsored  in  the  Department  of  Mechanical 
Engineering  and  Applied  Mechanics  at  The  University  of  Michigan 
b*  the  Army  Research  Office  during  the  period  1  January  1984 
through  31  August  1987. 

PROBLEM  STATEMENT 

During  this  program  of  research  we  have  explored  several 
topics  dealing  with  the  use  of  hybrid  optical/digital  methods  for 
recording,  visualizing  and  measuring  complex  flow  fields.  The 
ultimate  intent  of  this  research  was  to  develop  a  system  based  on 
digital  holographic  interferometry  and  computer  tomography  to 
determine  the  distribution  of  density  in  the  cross  section  of 
complex  three-dimensional  aerodynamic  flows.  We  have  been 
successful  in  this  endeavor. 

Our  work  had  three  basic  components.  The  first  was  the 
development  of  a  hybrid  system  to  record  interferometric  data  for 
subsequent  tomographic  analysis.  This  system  used  digital 
holographic  interferometry  (also  known  as  quasi-heterodyne 
holographic  interferometry  or  discrete  phase-step  holographic 
interferometry) .  The  system  had  to  be  capable  of  recording  the 
data  with  a  pulsed  laser,  had  to  record  data  for  a  large  number 
of  different  viewing  directions  simultaneously  and  had  to  be 
interpretable  in  a  relatively  automated  manner. 


The  second  component  of  our  work,  dealt  with  development  of 
reliable  methods  for  tomographic  reconstruction  of  flows  which 
give  rise  to  significant  refraction,  or  ray  bending,  of  the  laser 
light  used  to  probe  the  flow.  Ordinary  computer  tomography  is 
based  on  algorithms  which  assume  that  probing  rays  are  straight 
lines.  In  flows  containing  shocks  or  other  steep  gradients  this 
assumption  is  not  valid  because  the  consequent  variation  of 
refractive  index  causes  the  probing  rays  to  bend.  The  resulting 
problem  of  tomographic  reconstruction  is  both  nonlinear  and  ill- 
determined.  A  thorough  study  of  both  the  optical  and 
computational  aspects  of  this  problem  was  conducted  and  the 
techniques  developed  were  compared  in  detail  with  previously 
developed  methods  appropriately  modified  to  apply  to  aerodynamic 
measurements . 

The  third  component  of  our  research  was  the  establishment  of 
an  appropriate  computer  tomography  code  for  the  analysis  of 
complex  flow  data.  Although  the  literature  of  computer 
tomography  is  extensive,  exisiting  codes  did  not  appear  to  be 
satisfactory  for  our  purposes.  Therefore  we  set  out  to  develop  a 
new  technique  for  iterative  reconstruction  of  such  fields. 

Physical  experiments  were  planned  and  carried  out  using  a 
turbulent  free  jet  of  helium  as  the  test  flow.  Although  not 
originally  anticipated,  digital  holographic  interferometry  turned 
out  to  form  the  basis  of  a  useful  new  flow  visualization  system. 
This  avenue  of  research  was  also  explored  and  reported. 


In  this  section  we  outline  the  key  results  obtained  during 
the  research  program.  In  most  cases  the  results  have  been  fully 
documented  in  journal  articles,  which  are  referenced  in  the  text 
of  this  report. 

3 . 1  Digital  Holographic  Interferometry 

3.1.1  System  Description 

Digital  interferometry  is  a  recently  developed  hybrid 
optical-digital  metrology  technique  combining  two-exposure 
holographic  interferometry  with  digital  image  acquisition  and 
computer  processing  to  determine  the  interferometric  phase 
directly  from  a  set  of  image  irradiance  measurements.  The 
technique  is  similar  to  heterodyne  holographic  interferometry  in 
that  both  manipulate  the  interf erograms  phase  in  a  known  manner 
to  determine  its  magnitude.  Because  this  technique  requires  only 
the  recording  of  a  sequence  of  irradiance  values  of  each  pixel  in 
an  image  grid,  it  is  particularly  well  suited  to  recording  by 
vidicon  or  photodiode  arrays  and  processing  by  digital  computer. 
It  bypasses  many  of  the  problems  associated  with  analysis  of 
interferometric  fringe  patterns  and  readily  yields  the  sign  as 
well  as  magnitude  of  phase  shifts. 

In  this  type  of  interferometry,  a  two  exposure  holographic 
interf erogram  is  recorded.  One  exposure  is  made  of  the  flow  of 
interest  and  a  second  exposure  is  made  with  the  fluid  quiescent. 
Physically  different  reference  waves  are  used  for  each  of  the  two 
holographic  recordings.  At  the  time  of  reconstruction,  the  phase 


of  one  of  these  waves  is  changed  by  a  known  amount.  Such 


discrete  phase  changes  are  made  at  least  three  times  and  the 


consequent  change  of  irradiance  of  each  pixel  of  the  image  is 


recorded.  Elementary  mathematical  manipulation  enables  one  to 


calculate  the  unknown  phase  shift  (optical  pathlength  change) 


caused  by  the  flow  at  each  image  point.  More  than  three 


measurements  can  be  made  thereby  providing  computational 


redundancy  and  supressing  some  errors.  The  usual  complicated 


fringe  counting  procedure  is  replaced  by  a  simple,  computational 


sorting  operation. 


In  this  research  our  ultimate  objective  was  to  make 


tomographic  measurements  of  a  turbulent  jet;  therefore,  our 


system  needed  to  instantaneously  record  interferograms 


corresponding  to  several  different  viewing  directions.  Our 


system  was  based  on  the  use  of  a  pulsed  ruby  laser.  In  the 


course  of  the  experimentation  it  became  clear  that  this  laser 


needed  to  be  modified  in  order  to  perform  satisfactorily.  A 


temperature-stabilized  intracavity  etalon  was  installed  and 


thereafter  the  laser  produced  sufficient  coherence  and  mode 


structure  to  perform  the  experiments. 


The  test  section  was  a  rectangular  area  through  the  center 


of  which  a  free  helium  jet  flowed.  In  order  to  maximize  the 


number  of  viewing  angles  available,  two  sides  of  this  rectangular 


test  section  were  glassed  difusers  illuminated  from  behind  by  the 


pulsed  laser.  The  other  two  walls  of  the  rectangle  held  glass 


photographic  plates  on  which  the  holograms  were  recorded.  In 


order  to  get  these  holograms  as  close  as  possible  to  the  jet 


thereby  maximizing  the  solid  angle  of  view  available  one  of  these 


holograms  was  a  traditional  transmission  hologram  while  the  other 
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was  a  reflection  hologram.  A  reflection  hologram  is  one  in  which 
the  reference  wave  and  object  wave  enter  opposite  sides  of  the 
holographic  plate  and  the  holographically  reconstructed  object 
wave  is  formed  by  reflection  from  the  developed  plate  rather  than 
by  transmission.  The  development  of  this  system  required  a 
significant  amount  of  experimentation  with  reflection  holograms. 

A  variety  of  chemical  processing  techniques  for  development  and 
bleaching  of  reflection  holograms  was  tried  and  comparatively 
evaluated . 

A  reconstruction  system  using  a  helium-neon  laser  was 
constructed  on  a  separate  optical  table.  This  system  included  a 
mirror  mounted  on  a  piezo-electric  crystal  with  a  feedback 
control  system  for  varying  the  phase  by  known  amounts  and  holding 
that  phase  stable  during  the  reconstruction  process.  Irradiance 
data  were  recorded  using  a  high  quality  Cohu  video  camera. 
Digitized  images  recorded  in  this  manner  were  stored  in  a  local 
LSI  11/23  microcomputer  which  was  also  used  to  analyze  the  phase 
shifts  and  to  display  computed  images.  Tomographic 
reconstructions  were  carried  out  either  on  The  University  of 
Michigan's  mainframe  IBM  3090  computer  or  on  an  Apollo  computer 
workstation . 

This  basic  system  and  technique  is  described  in  more  detail 
in  a  paper  entitled  "Digital  Interferometry  for  Flow 
Visualization"  which  is  included  in  the  appendix  of  this  report 
and  will  be  more  fully  described  in  a  journal  article  currently 
in  preparation. 


1 


y  \*v 

t  ■  -  n  »1  1a  M.  ^  V.  aA»A  1  -  W_  A  V_  V-  V*_  V 


Flow  Visualization 


An  interesting  and  very  useful  sidelight  of  the  research 
program  was  the  development  of  a  flow  visualization  technigue 
based  on  digital  holographic  interferometry.  Recall  that  the 
optical  pathlength  change,  or  phase  shift,  at  each  point  in  the 
image  of  the  flow  is  a  computed  value.  This  is  different  than 
ordinary  interferometry  which  directly  results  in  an  analog 
interference  fringe  pattern.  It  therefore  is  possible  to  display 
results  in  a  number  of  ways.  Ordinarily,  one  would  display 
contours  of  constant  phase  shift,  that  is  calculate  an 
interf erogram .  However,  if  one  displays  a  simple  gray  scale 
image,  for  example  setting  a  value  of  zero  to  the  largest  phase 
shift  occurring  and  a  value  of  one  to  an  undisturbed  region  of 
the  flow,  a  remarkably  clear  visualization  of  the  flow  can  be 
displayed.  Such  images  are  very  similar  in  appearance  to  images 
formed  by  techniques  such  as  laser  induced  fluorescence.  The 
properties  of  such  images  are  similar  to  traditional  shadowgraphs 
except  they  are  sharply  in  focus  in  the  object  plane  of  the 
camera  used  to  record  the  image.  Such  images  may  be  particularly 
useful  for  relatively  weak  flows,  because  the  digital 
interferometric  system  is  very  accurate  and  can  be  used  to  record 
images  that  would  result  in  less  than  a  single  fringe  in  an 
ordinary  interferogram .  This  technique  is  recorded  briefly  in  a 
paper  in  the  appendix  to  this  report  as  well  as  in  the  PhD  thesis 
of  David  W.  Watt,  and  will  also  appear  in  a  journal  article 
currently  under  preparation. 


3.2  Tomographic  Code:  Consistent  Iterative  Convolution 


In  the  current  study  we  attempted  to  compute  cross  sectional 
images  of  a  turbulent  helium  jet.  This  flow  is  relatively 
complex  and  requires  a  considerably  greater  amount  of  data  than 
have  been  used  for  previous  reconstructions  in  this  laboratory. 
The  question  of  how  best  to  effect  the  tomographic  reconstruction 
is  not  a  simple  one  and  was  investigated  in  some  detail  during 
this  research.  We  settled  on  the  use  of  iterative  convolution. 
This  means  that  the  fundamental  reconstruction  technique  is  the 
classical  convolution  algorithm.  However,  this  algorithm  assumes 
that  experimental  data  are  available  for  a  full  set  of  equally 
spaced  viewing  directions  distributed  over  180  degrees.  This 
requirement  could  not  be  met  in  the  experimental  set  up  developed 
for  this  research,  nor  is  it  likely  to  be  in  most  fluid  dynamic 
experiments.  Hence  we  developed  an  iterative  technique  in  which 
the  missing  data  are  effectively  replaced  computationally.  Such 
an  iterative  convolution  technique  was  originally  proposed  by 
this  laboratory  under  a  previous  ARO  contract.  However,  the 
basic  iterative  convolution  technique  was  not  sufficiently 
accurate  to  handle  the  complexity  occuring  in  this  experiment, 
and  its  convergence  properties  were  poor  and  poorly  understood. 

In  order  to  analyze  the  data  from  this  experiment  we 
extended  our  previous  work  by  developing  a  technique  of 
constrained  iteration  between  the  estimated  image  of  the  flow 
and  its  Radon  transform.  (The  Radon  transform  is  the  appropriate 
mathematical  term  for  the  set  of  line  integrals  of  flow  density 
obtained  interf erometrical ly . )  We  found  that  for  certain  test 
functions  used  to  study  this  technique,  the  estimated  image 


diverges  from  the  actual  function.  This  divergence  was  shown  to 
result  from  the  algorithm's  failure  to  enforce  consistency  between 
the  image  estimate  and  its  measured  Radon  transform.  Such 
consistency  in  the  new  approach  is  enforced  using  a  routine  based 
on  the  direct  inversion  formula  and  on  decomposing  the  image  into 
its  measured  and  missing  components.  This  routine  was  used  to 
develop  a  new  iterative  algorithm  which  converges  absolutely  for 
all  test  functions  studied. 

The  consistent  iterative  convolution  technique  and  its 
behavior  are  described  in  detail  in  a  journal  article  currently 
under  preparation  and  also  in  the  doctoral  thesis  of  David  W. 

Watt . 

3 . 3  Measurement  of  Turbulent  Helium  Jets 

The  total  system  for  holographic  recording,  digital  image 
aquisition  and  processing,  and  tomographic  reconstruction  was 
developed  and  studied  by  examining  free  helium  jets  in  air.  The 
jets  studied  had  Froude  numbers  varying  from  160  to  37,500  and 
Reynolds  numbers  varying  from  250  to  1250. 

Prior  to  conducting  tomographic  reconstructions,  we  examined 
flow  visualization  images  of  several  of  the  jets.  These  images 
tended  to  confirm  contentions  previously  appearing  in  literature 
that  such  turbulent  jets  are  organized  into  large  scale  motions 
whose  length  scale  is  of  the  order  of  the  local  flow  width.  The 
tomographic  images  of  the  cross  sections  of  the  jets  clearly 
confirm  the  presence  of  nearly  unmixed  ambient  fluid  within  the 
jet  boundaries  and  indicated  that  these  inclusions  are 
responsible  for  steep  concentration  gradients  within  the  jet. 


These  images  also  seem  to  show  the  presence  of  vortical  motions 
in  the  off-center  jet  regions  which  result  in  medium  scale  (half 
the  jet  width  or  less)  inclusions  of  ambient  fluid.  They  also 
seem  to  confirm  that  the  jet  scalar  field  is  divided  into  regions 
of  largely  constant  concentration.  It  must  be  recognized  that  we 
did  not  conduct  a  detailed  study  of  this  turbulent  flow  and  that 
the  system  we  have  constructed  can  "freeze"  the  entire  three- 
dimensional  flow  field  at  an  instant  of  time,  but  does  not 
document  the  time  dependent  behavior  of  the  jet.  The  images  also 
are  imperfect  and  contain  certain  experimental  and 
reconstructional  artifacts.  Some  of  these  difficulties  stem  from 
the  fact  that  we  used  lasers  of  different  frequencies  for 
recording  and  reconstruction.  Although  we  designed  the  setup  to 
minimize  apparent  deformations  caused  by  this  wavelength  shift, 
they  are  not  eliminated  entirely.  Nonetheless,  we  consider  these 
experiments  to  be  highly  successful  and  to  show  the  practicality 
of  making  instantaneously  three-dimensional  measurements  in 
complex  flows. 

3 . 4  Tomography  of  Strongly  Refracting  Flows 

When  strong  density  gradients  occur  within  the  flow  field 
under  examination,  the  probing  optical  rays  may  be  significantly 
bent.  In  this  case  the  well-known  computer  tomography 
algorithms,  which  assume  probing  rays  to  be  straight  lines,  are 
not  applicable.  During  this  research  program,  we  carried  out 
an  extensive  study  of  tomographic  reconstruction  in  the  presence 
of  such  strong  refraction.  Related  problems  are  somewhat  better 
explored  in  the  context  of  acoustical  and  ultrasonic  problems. 


Distinct  differences  arise  when  optical  problems  are  considered. 
In  this  research  we  posed  the  problem  in  the  optical  context. 

The  resulting  mathematical  problem  is  highly  nonlinear,  and  we 
approached  it  by  two  iterative  techniques  and  a  perturbation 
analysis . 

The  first  procedure  is  called  straight  line  inversion 
with  modified  data  (SLIM)  and  is  based  on  previous  work  by  Bates 
and  McKinnon  and  by  Cha  and  Vest.  This  is  essentially  a  method 
of  successive  approximations  in  which  an  initial  estimate  of  the 
field  is  computed  by  ignoring  the  refraction  effects  and  then 
iteratively  corrected  by  using  computational  ray  tracing  to 
effectively  calculate  an  interf erogram  which  can  be  compared  with 
the  measured  data.  The  second  iterative  technique  is  refered  to 
as  the  curved  ray  algebraic  inversion  (CRAI) .  This  is  based  on 
earlier  work  of  Johnson  and  of  Schomberg.  In  this  procedure  the 
current  estimate  of  the  index  of  refraction  is  used  to  generate  a 
system  of  algebraic  equations  where  the  unknowns  are  corrections 
to  this  estimate  at  each  pixel.  In  a  sense  it  is  a  modified 
version  of  the  well-known  ART  procedure. 

The  third  procedure  studied  was  a  perturbation  approach 
which  does  not  involve  iteration.  Such  an  approach  was  orginally 
developed  for  ultrasonic  problems  by  Norton  and  Linzer.  It 
assumes  that  the  index  of  refraction  at  each  point  varies  only 
slightly  from  its  value  in  the  surroundings  and  that  the  ray 
tradjectory  deviates  very  little  from  a  straight  line.  We 
modified  the  original  analysis  to  properly  represent  data 
acquisition  by  holographic  interferometry  including  an  imaging 
system.  This  results  in  a  very  simple  computation.  We  found 


that  both  the  the  SLIM  and  CRAI  techniques  tend  to  initially 
converge  to  a  more  accurate  estimated  reconstruction  and  then 
diverge.  Through  our  study  we  gained  some  understanding  of  the 
causes  of  this  divergence.  We  also  found  that  although  the  CRAI 
approach  is  more  difficult  to  implement  and  slower  than  SLIM,  it 
is  to  be  preferred  in  that  it  generally  produces  less  pronouned 
divergence . 

The  primary  conclusion  drawn  from  this  work,  however,  is 
that  under  most  circumstances  the  simple  perturbation  technique 
performs  admirably.  Particularly  in  the  case  of  aerodynamic 
flows  where  strong  divergence  from  straight  optical  paths  is 
usually  small,  this  would  appear  to  be  the  method  of  choice. 
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APPENDIX 


This  appendix  contains  copies  of  journal  articles  based  on 
this  research  which  have  been  published  to  date.  Two  more 
comprehensive  articles  are  in  preparation. 


I 


► 


vi 


15 


» 

y 


T’jf^’VO^ 


Reprinted  from  Applied  Optica.  Vol.  2/,  page  1089,  December  I.  1985 
Copyright  ©  1985  by  the  Optical  Society  of  America  and  reprinted  by  permission  of  the  copyright  owner. 

Tomography  for  properties  of  materials  that  bend  rays: 
a  tutorial 


Charles  M.  Vest 


When  tomography  is  performed  with  electromagneticor  acoustical  radiation,  refraction  may  cause  sufficient 
bending  of  the  probing  rnva  that  ordinary  reconstruction  algorithms,  which  are  based  on  the  assumption  of 
straight  rays,  do  not  yield  accurate  results  The  resulting  problem  of  reconstructing  the  refractive-index 
distribution  of  an  object  from  time  of  flight  or  optical  path  length  data  is  nonlinear.  Various  approaches  to 
solving  this  problem  approximately  have  been  proposed  and  subjected  to  modest  numerical  studies.  These 
include  iterative  algorithms  and  techniques  based  on  linearized  inverse  scattering  theory.  One  exception  is 
the  case  of  axisymmetric  objects  for  which  an  exact  solution  is  known. 


I.  Introduction 

The  classical  technique  of  computed  tomography  is 
based  on  the  assumption  that  the  probing  or  emitted 
rays  of  radiation  are  straight  lines.  This  assumption 
leads  to  the  now  well-known  problem  of  reconstruction 
of  some  density  function  from  experimentally  collect¬ 
ed  values  of  its  line  integrals.  This  problem  is  linear 
and  can  be  approached  with  the  mathematics  of  the 
Radon  transform  and  associated  techniques.  Here  we 
are  concerned  with  cases  where  the  probing  or  emitted 
rays  are  bent  bv  refraction  associated  with  gradients  of 
speed  (refractive  index)  within  the  object.  Hence  the 
density  function  must  be  reconstructed  from  mea¬ 
sured  values  of  path  integrals  along  generally  unknown 
curved  paths.  This  problem  is  inherently  nonlinear 
and  is  not  generally  associated  with  known  mathemati¬ 
cal  transforms. 

In  this  brief  tutorial  paper,  I  will  discuss  several 
aspects  of  the  problem  of  computed  tomography  of 
objects  that  bend  rays  and  will  outline  some  of  the 
approaches  which  have  been  taken  to  reconstruct  such 
objects.  These  problems  are  associated  primarily  with 
ultrasonic  and  optical  tomography.  In  most  applica¬ 
tions  of  tomography  using  x  rays  and  nuclear  sources, 
the  wavelengths  are  too  small  for  such  phenomena  to 
be  important.  My  use  of  the  term  rays  in  the  title  and 
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Sec.  I  displays  a  personal  tendency  to  think  about  the 
problem  in  terms  of  ray  optics.  Other  approaches 
summarized  herein  are  based  on  wave  theory  and  in¬ 
verse  scattering. 

In  applications  to  both  ultrasonics  and  optics  (or 
other  domains  of  electromagnetics),  the  objective  usu¬ 
ally  is  to  reconstruct  either  the  distribution  of  an  at¬ 
tenuation  (or  emission)  coefficient  or  of  refractive  in¬ 
dex.  In  the  case  of  ultrasonics  the  measured  data  are 
time  of  flight,  intensity,  or  complex  amplitude  (real 
amplitude  and  phase).  In  the  case  of  optics  the  mea¬ 
sured  data  are  complex  amplitude  or  intensity.  The 
intensity  may  in  fact  be  a  fringe  pattern  when  the  data 
are  gathered  interferometrically. 

Particularly  in  the  case  of  ultrasonics,  refraction 
may  cause  secondary  problems — if  one  is  measuring 
attenuation,  refraction  can  cause  divergence  of  the 
beam  as  it  travels  from  the  transmitter  to  the  receiver, 
thereby  giving  rise  to  an  apparent  attenuation.  Such 
problems  are  compounded  if  the  signal  detection  is 
phase  sensitive,  because  phase  cancellation  effects  can 
occur.  Problems  of  this  class  have  been  studied  by 
several  authors,  and  a  variety  of  approximate  correc¬ 
tion  schemes  have  been  developed;  see,  for  example, 
Farrell,1  Pan  and  Liu,2  Klepper  et  al and  Itoh  et  al.* 

Refractive  effects  can  be  distributed  continuously 
throughout  the  object  volume,  or  they  can  be  discrete, 
as  at  an  object  boundary.  For  example.  Eberhard5  has 
dealt  with  ultrasonic  time-of-flight  tomography  for 
nondestructive  testing  of  turbine  blades.  To  over¬ 
come  errors  due  to  very  strong  bending  of  the  rays  at 
the  object  boundary,  he  encapsulated  the  blade  in  a 
rectangular  block  of  material  whose  refractive  index 
nearly  matched  that  of  the  blade.  Because  of  this 
simple  rectangular  cross  section,  an  algorithm  could 

1  December  1985  /  Vol  24.  No  23  /  APPLIED  OPTICS  4089 


. V.V.V.VV 


easily  he  developed  to  correct  for  the  effects  of  refrac¬ 
tion  at  its  boundary. 

In  this  paper,  we  are  primarily  concerned  with  re¬ 
fractive  effects  that  are  distributed  continuously 
throughout  the  object.  Data  are  assumed  to  be  either 
time  of  flight  or  optical  path  length. 

II.  Axisymmetric  Objects:  Reconstruction  from  Path 
Length  Data 

It  is  useful  to  consider  tomography  based  on  optical 
path  length  measurements  of  axisymmetric  objects. 
This  case  has  an  analytical  solution  under  modestly 
constrained  conditions.  It,  therefore,  gives  insight 
into  the  more  general  problem.  Also,  experience  with 
interferometry  of  axisymmetric  objects  shows  the  im¬ 
portant  role  played  bv  the  physics  of  the  data  gather¬ 
ing  method,  imaging  in  this  case. 

A  typical  ray  traversing  a  cylindrical  object  with 
axisymmetric  refractive  index  distribution  n(r)  is 
shown  in  Fig.  1.  If  refractive  effects  were  negligible, 
the  rays  through  this  object  would  be  straight  lines. 
We  then  would  measure  the  optical  path  length  (line 
integral)  of  all  rays  in  a  given  direction,  say  parallel  to 
the  v  axis.  This  path  length,  which  we  denote  as  <p{x), 
is  given  by 


(t(X) 


f°  f[r)rdr 
,  (r2  -  x2)1'2  ’ 


(1) 


which  is  an  Abel  integral  of  n(r).  The  inversion  of  Eq. 
(1)  is  well  known: 


Hr)  - 


I  l'°  ld<t>/dx)dx 
lx2”—  r2)^2 


(2) 


Thus  by  measuring  one  projection  we  can  reconstruct 
the  refractive-index  distribution.  Note  that  this  is  a 
linear  problem.  A  large  number  of  numerical  algo¬ 
rithms  have  been  devised  to  approximate  this  recon¬ 
struction  when  data  are  discrete. 

Now  let  us  return  to  the  case  where  refraction  is 
strong.  By  definition,  the  optical  path  length  of  the 
ray  shown  in  Fig.  1  is 


>t>  =  |  tuis. 


where  ds  is  the  differential  length  of  the  ray.  Al¬ 
though  we  do  not  know  the  rav  curve,  we  do  know  that 
it  is  governed  by  the  rav  equation 


which  has  a  simple  solution,  for  the  axisymmetric  case, 
termed  Bouguer's  formula": 

rr?< r  I  sind  I  =  p  (M 

Note  that  the  rays  are  parametrized  by  the  constant,  p. 
sometimes  called  the  impact  parameter.  By  introduc¬ 
ing  a  new  variable. 

•)  —  r  f  H  r  |  (8  I 

together  with  geometric  relations  derived  from  Fig  1. 
the  following  expression  for  the  optical  path  length  can 
be  found: 

■:'r>  -  I  ■ 
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Fig.  1.  Optical  ray  passing  through  a  strongly  refracting  axisym¬ 
metric  object. 


It  is  quite  interesting  to  note  that  although  the  prob¬ 
lem  we  are  solving  is  highly  nonlinear,  Eq.  (7)  is  the 
Abel  integral  of  the  quantity  [rjcf  ln(r)/dr/];  hence  its 
solution  is  of  the  form  of  Eq.  (2): 


d  Inr\  _  1  I*  ( d<t>/dp\dp 

dV)  7  I,  (p2  -  ,2)>'2  ‘ 


18) 


If  desired,  this  can  be  integrated  to  give 


[ip  cosh-  { 

p\d<b  dp 

M.  v 

w)  dp  p 

Thus  given  the  set  of  path  integrals  and  the  values  of 
the  impact  parameter  (i.e.,  the  entrance  or  exit  direc¬ 
tion  and  the  refractive  index  at  the  boundary),  the 
refractive-index  distribution  is  reconstructed  implicit¬ 
ly  and  can  be  determined  as  long  as  the  quantity  rn(r) 
is  a  monotonic  function  of  r.  Further  discussion  of  this 
type  of  reconstruction  can  be  found  in  Bullen,7  Phin- 
ney  and  Anderson,8  and  Vest9  in  the  context  of  geo¬ 
physical  acoustics,  radio  exploration  of  planetary  at¬ 
mospheres,  and  interferometry  of  gases,  respectively. 

An  interesting  effect  occurs  when  axisymmetric 
strongly  refracting  objects  are  studied  by  interferome¬ 
try.  In  this  case  the  data  are  recorded  in  the  form  of  a 
fringe  pattern  formed  by  the  interference  of  a  refer¬ 
ence  plane  wave  with  an  initially  plane  wave  that  has 
been  distorted  by  travelling  through  the  object.  From 
this  fringe  pattern  one  can  determine  the  distribution 
of  phase  across  the  test  wave.  If  the  refraction  is 
negligible,  all  rays  involved  are  essentially  straight  and 
parallel  to  each  other.  We  again  have  a  simple  Abel 
problem  which  can  be  solved  for  j(r)  =  n(r)  —  rip,  where 
n,,  is  the  uniform  refractive  index  of  the  medium  sur¬ 
rounding  the  object. 

Now  suppose  that  the  object  strongly  refracts  the 
test  wave,  as  in  Fig.  1.  It  then  is  important  to  intro¬ 
duce  an  imaging  lens  between  the  object  and  the  plane 
in  which  the  fringe  pattern  is  to  be  observed.  This  lens 
can  be  focused  on  anv  plane  parallel  to  the  observation 
plane  It  has  been  shown  by  extensive  numerical  ex¬ 
periments, 1  that  if  the  lens  is  focused  on  the  center 
plane  of  the  object,  and  if  one  applies  the  Abel  inver¬ 
sion  to  the  data  as  if  no  refraction  had  occurred,  the 
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reconstruction  will  be  very  accurate,  even  if  refraction 
is  extremely  strong.  If  another  plane,  away  from  the 
center,  is  imaged,  this  procedure  will  not  give  a  good 
result.  Hence  imaging,  or  an  equivalent  data  process¬ 
ing  operation,  essentially  eliminates  the  problems  as¬ 
sociated  with  strong  refraction  in  interferometry  of 
axisymmetric  objects. 

III.  Asymmetric  Objects:  Reconstruction  from  Path 
Length  Data 

Consider  the  data  gathering  system  shown  in  Fig.  2. 
It  is  assumed  that  the  field  of  interest  is  two-dimen¬ 
sional  and  lies  in  a  circular  region.  Data  are  collected 
by  interferometry;  the  fringe  pattern  is  formed  in  the 
image  plane  by  an  imaging  system  focused  on  a  plane 
located  a  distance  rf  from  the  center  of  the  circle.  The 
unknown  refractive-index  field  is  represented  by 
n(r,</>),  the  optical  axis  is  aligned  so  that  in  the  absence 
of  refraction  the  projection  direction  would  be  speci¬ 
fied  by  the  angle  0.  and  the  optical  path  length  value  is 
recorded  at  a  point  P  in  the  image  plane,  which  is  the 
conjugate  of  the  point  Pin  the  object  space.  As  shown 
in  the  figure,  two  rays  meet  and  interfere  at  point  P: 
one  travels  a  curved  path  through  the  object  because  of 
refraction  due  to  nlr,0),  and  the  other  is  a  reference  ray 
which  travels  along  a  straight  path  through  a  medium 
with  uniform  refractive  index  n,). 

Using  elementary  geometric  optics  it  is  seen  that  the 
corresponding  optical  path  length  difference  is  given 
by 

<K  --  __  __ 

A</>(p.H)  *  J  ri(r,it>)ds  +  fi(1(/JC  -  I)E  -  F.F) 

=  n„|.  •  (10) 

Cha  and  Vest11  have  termed  this  the  optical  path 
length  transform.  INote  that  if  ray  bending  is  negligi¬ 
ble.  it  reduces  to  the  line  integral  transform,  A<j>  = 
P\nir,<t>);  n,,|  to  which  the  usual  techniques  of  comput¬ 
ed  tomography  can  be  applied.!  Two  approaches  have 
been  proposed  to  obtain  reconstructions  from  data  of 
this  general  type:  iterative  algorithms  and  perturba¬ 
tion  analysis.  Both  take  as  a  starting  point  an  initial 
reconstruction  made  hy  ignoring  refraction.  To  dis¬ 
cuss  these,  a  deviation  function  which  is  the  difference 
between  the  path  length  transform  and  the  line  inte¬ 
gral  transform  of  the  object,  or  of  an  estimate  of  the 
object,  is  defined: 

l Up, (I (  =  Ailp.tl)  -  A-/>ip.«l  Ul) 

The  iterative  procedure  was  devised  bv  Cha' and  is 
discussed  in  detail  in  Ref.  11.  The  algorithm  is  as 
follows: 

(i)  Make  an  initial  estimate  of  the  deviation  func¬ 
tion  Dip  ft),  where  i  =  0.  (Note  that  we  need  not  start 
with  D, i  =  0;  if  one  has  some  a  priori  knowledge  about 
the  general  structure  of  the  object ,  convergence  mav  he 
speeded  hy  guessing  the  rough  structure  of  the  devi¬ 
ation  function.) 

( ii)  Calculate  the  corresponding  estimate  of  the  line 
integral  translorm: 


Fig.  2.  Formation  of  an  interferogram  of  an  asymmetric  refracting 
object. 


A«,(p.fl)  -  Aiip.6)  -  D,(p.»).  (12) 

(iii)  Reconstruct  the  object  approximately  by  com¬ 
puting  the  inverse  line  integral  transform; 

n,(r,0)  -  n„  =  P~'(Ai>,).  (13) 

(iv)  Using  computational  ray  tracing,"'  calculate 
the  path  length  transform  of  the  estimated  distribu¬ 
tion: 

A<A,(p.»)  =  /’In.fr,^);  n0|.  (14) 

(v)  Compute  a  new  estimate  of  the  deviation  func- 


(vi)  Return  to  step  (ii)  and  continue  the  process 
iteratively  until  some  measure  of  D,(p,<p)  is  sufficiently 
small. 

Cha  and  Vest1 1  have  applied  this  algorithm  to  sever¬ 
al  numerical  experiments  in  which  data  computed  for  a 
specified  object  are  used  as  input  to  the  algorithm,  and 
its  convergence  toward  an  accurate  reconstruction  was 
studied.  The  algorithm  did  indeed  produce  rather 
accurate  reconstructions  in  the  specific  cases  studied; 
however,  some  operator  interaction  was  needed  to  de¬ 
tect  computational  ray  crossing.  Path  length  data 
contaminated  by  this  effect  were  eliminated.  They 
also  applied  a  modified  form  of  this  algorithm  to  ex¬ 
perimental  data  obtained  from  interferometry  from 
strongly  refracting  electrochemical  boundary  layers. 
Related  work  on  iterative  correction  for  ray  bending 
has  been  carried  out  by  Greenleaf  and  Johnson  (see, 
e.g.,  Ref.  14)  and  Glover  and  Sharp.15 

Another  iterative  approach  was  developed  by 
Schomberg.lfi  He  devised  a  modified  ART  procedure 
to  determine  refractive  index  at  discrete  pixels.  For 
each  iteration,  the  coefficients  of  the  algebraic  equa¬ 
tions  are  determined  by  ray  tracing  using  the  values  of 
refractive  index  from  the  previous  iteration.  This  pro¬ 
cess  is  continued  until  the  computed  and  measured 
path  integrals  are  in  good  agreement. 

McKinnon  and  Bates17  have  presented  another  iter¬ 
ative  algorithm  that  does  not  require  ray  tracing  at 
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each  iteration.  They  accomplish  this  by  the  less  time- 
consuming  task  of  solving  an  approximate  form  of  the 
eikonal  equation. 

To  summarize,  several  iterative  algorithms  have 
been  proposed.  Each  has  been  demonstrated  to  pro¬ 
duce  improved  reconstructions  in  a  small  number  of 
experiments  or  numerical  simulations.  On  the  other 
hand,  each  of  these  is  a  somewhat  ad  hoc  procedure  for 
which  no  proof  or  convincing  argument  for  conver¬ 
gence  under  general  conditions  has  been  given.  In 
fact,  McKinnon  and  Bates17  indicate  that  in  the  com¬ 
mon  format  for  recording  ultrasonic  data,  refracting 
objects  have  forbidden  regions  which  cannot  be  recon¬ 
structed,  and  that  it  may  seldom  be  feasible  to  improve 
significantly  images  beyond  those  reconstructed  bv 
ignoring  refraction.  Although  this  may  be  overly  pes¬ 
simistic  in  view  of  some  of  the  successful  applications 
to  smoothly  varying  objects,  it  is  clear  that  such  proce¬ 
dures  must  be  applied  with  caution,  and  that  more 
study  of  convergence  and  algorithm  behavior  is  need¬ 
ed. 

Although  iterative  procedures  are  an  obvious  ap¬ 
proach  to  the  problem  under  consideration,  thev  are 
inherently  slow  and  expensive  because  ray  tracing  is 
computationally  time-consuming.  Hence  there  is  mo¬ 
tivation  to  seek  more  direct  approaches.  Norton  and 
Linzer18  have  carried  out  an  extensive  study  of  the 
application  of  perturbation  analysis  to  the  problem  of 
reconstruction  of  objects  that  bend  rays.  In  essence, 
they  also  start  with  a  deviation  function  representing 
the  difference  between  data  received  from  the  refract¬ 
ing  medium  and  that  which  would  be  received  if  the 
probing  rays  were  straight  lines.  However,  they  then 
seek  a  single  operation  which  can  be  directly  applied  to 
this  deviation  function  to  generate  an  accurate  approx¬ 
imate  reconstruction. 

Because  Norton  and  Linzer's  analysis  is  designed  to 
be  applied  directly  to  time-of-flight  ultrasonic  data 
with  no  imaging  or  interferometry,  it  is  convenient  to 
use  their  notation,  which  differs  from  that  above.  It  is 
assumed  that  the  ultrasonic  refractive  index  deviates 
from  unity  by  only  a  small  amount: 

'Hr)  —  1  +  Wi(rl.  (16) 

where  t«  1  is  a  small  parameter.  Correspondingly,  as 
indicated  in  Fig.  3,  the  actual  ray  path  £  deviates 
slightly  from  the  straight  path  L.  If  R  is  a  vector 
extending  from  the  transmitter  T  to  the  receiver  R,  the 
time  of  flight  through  the  object  is 


r,iRl=  rilrlds.  (IT) 

!/  'Ri 

where  the  subscript  reminds  one  that  this  time  is  asso¬ 
ciated  with  the  curved  path  £,  as  opposed  to  the 
straight  path  /..  Now  let  h(R)  be  an  estimate  of  the 
refractive  index  obtained  from  a  straight  line  algo¬ 
rithm.  i.e.,  a  reconstruction  obtained  by  neglecting 
refraction, 

T,  ■  K  •  =  I  Oi  Rlr/s.  (IS) 

/!  ■  «' 

A  time -delay  correction,  or  deviation  function,  is  de¬ 
fined  as 
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Pig.  3.  Rav  passing  through  a  mildly  refracting  object.  It  origi¬ 
nates  at  transmitter  T  and  impinges  on  receiver  R. 


e(R,A)  =  7^  (R)  -  T,(R).  (19) 

The  objective  of  this  technique  is  to  derive  a  directly 
calculable  expression  for  e(R,n). 

Because  the  ray  path  is  assumed  to  deviate  only 
slightly  from  the  straight  line  L,  and  because  the  re¬ 
fractive  index  along  £  differs  only  slightly  from  that 
along  the  line,  each  can  be  expanded  in  power  series  in 
the  perturbation  parameter  t  as 

y(x>  -  r/(x>  +  0(<2)  (20) 


n(x.y)  =  1  +  (hn(x)  +  rf(x)h^x)  +  .  .  . ,  (21) 

where  y(.t)  is  the  trajectory  of  the  ray  passing  through 
O  and  R  and  where 

hn(x)  =  h(x.O).  h,(i)  =  (<Vj/<5y)|,.fl.  (22) 

By  applying  Fermat's  principle  and  the  calculus  of 
variations.  Norton  and  Linzer18  minimize  Ti  with  re¬ 
spect  to  y(x )  and  obtain  the  explicit  expression  for  the 
deviation  function  e: 


where  f(x)  is  the  solution  of 


subject  to 


d2(/dx7  =  hjx) 


/(())  -/(!)-«. 


There  are  several  assumptions  which  must  be  con¬ 
sidered  if  one  wishes  to  apply  this  technique.  First 
and  foremost,  because  it  is  based  on  a  perturbation 
analysis,  variations  of  refractive  index  must  be  small. 
Numerical  simulations  suggest  that  velocity  variations 
should  be  limited  to  about  ±5%.  which  suggests,  for 
example,  that  it  should  be  applicable  to  soft  tissue 
measurements.  Variations  of  refractive  index  also 
should  be  smooth.  The  application  of  ray  optics  must 
be  appropriate,  and  the  transducer  width  w  should  be 
sufficiently  small: 


O  V.- 


Vu7  «  a.  (26) 

where  a  is  the  characteristic  scale  of  the  inhomogenei¬ 
ties  of  refractive  index,  and  1  is  the  propagation  dis¬ 
tance  across  the  object. 

IV.  Methods  Based  on  Diffraction  and  Inverse  Scattering 
Theory 

Rather  than  utilizing  a  ray-optical  approach,  it  is 
possible  to  view  reconstruction  of  refractive  index  as  a 
problem  in  inverse  scattering  theory.  If  diffraction  is 
mild,  this  can  be  approached  by  using  the  Rvtov  or 
Born  approximations  to  the  wave  equation.  Iwata  and 
Nagata19  first  applied  this  approach  to  problems  of 
reconstruction  from  optical  data  in  1975.  Kenue  and 
Greenleaf20  recently  considered  the  use  of  the  Rytov 
approximation  for  reconstruction  from  ultrasonic 
time  -of-flight  data.  Their  approach  is  briefly  outlined 
below. 

Figure  i  depicts  a  plane  wave  of  radiation  impinging 
on  the  object  /(.t,_v)  of  interest  and  being  refractively 
distorted  as  it  propagates  through  it.  At  some  plane 
beyond  the  object  the  complex  amplitude 

('(r)  =  exp|iKn<Mr)|  (27) 

of  the  emerging  wave  is  measured  in  an  observation 
plane.  This  is  repeated  for  a  number  of  different 
propagation  directions.  It  is  assumed  that  the  Helm¬ 
holtz  wave  equation  applies: 

V2(’(r)  +  |Knn(r)|2(’(r)  =  0.  (28) 

The  eikonal  then  must  obey  the  relation 

iKn-'r2*  -  (V</>)2  +  n2  =  0.  (29) 

Remembering  that  the  variation  or  refractive  index 
nfr)  is  small,  it  is  appropriate  to  write 

n  -  1  +  n  | .  tf>  —  0ft  +  0,,  (30) 

where  <Pt\  is  the  eikonal  in  the  undisturbed  medium. 
For  the  entering  plane  wave, 

•/>„  =  K  ■  R.  CM) 

Integrating  across  the  object  we  obtain  the  line  integral 


’In— ol 

-  '/>  =  n^x.y)ds. 


When  Eq.  (30)  is  substituted  into  Eq.  (29),  we  obtain  to 
first  order 

^“0,  +  K„J0,  =  (33) 

For  a  lixed  frequence  K, >  and  direction  ».  the  2-D 
Fourier  transform  solution  ot  this  equation  is 

nU’.I’I  =  in  pxpjfi  A‘„  -  u  (.Y(l|0,(.Y„.r,OI.  <341 

where  u  and  r  are  the  transform  variables  in  the  direc¬ 
tion  of  K  and  normal  to  K.  respectively. 

The  complex  log  of  the  measured  values  of  t’(r)  is 
calculated  to  give  cn.  whose  I  I)  Fourier  transform  is 
then  computed  and  substituted  info  Eq  (31).  This 
gives  the  2-1)  transform  along  a  circle  in  the  2-1)  trans¬ 
form  plane.  Note  that  to  obtain  values  around  the 


U(x0.y) 


Fig.  4.  Complex  amplitude  of  a  wave  (initially  plane)  that  has 
propagated  through  a  refracting  object. 

entire  circle  requires  both  transmission  and  reflection 
data.  If  this  is  repeated  for  several  different  propaga¬ 
tion  directions,  the  entire  Fourier  transform  plane  can 
be  filled  in.  The  refractive-index  distribution  can 
then  he  computed  by  inverse  transformation. 

It  is  interesting  to  note  the  analogy  between  the 
above  result  and  the  central  section  theorem  of  ordi¬ 
nary  computed  tomography. 

Recent  work  by  Devaney2122  has  treated  this  type  of 
inverse  problem  for  weakly  scattering  objects  using  the 
Born  approximation.  Again  the  fundamental  result  is 
that  by  computing  1-D  Fourier  transforms  of  the  com¬ 
plex  amplitude  of  the  scattered  field  it  is  possible  to 
determine  the  2-D  transform  of  the  object’s  refractive 
index  along  circular  arcs.  He  has  generated  a  filtered 
backpropagation  algorithm  that  operates  in  a  manner 
somewhat  analogous  to  the  well-known  filtered  back- 
projection  algorithm  of  ordinary  computed  tomogra¬ 
phy.  This  interesting  and  important  techn  que  is  dis¬ 
cussed  in  a  paper  by  Devaney.23 

V.  Closure 

In  this  paper  I  have  outlined  some  of  the  techniques 
which  have  been  proposed  for  reconstruction  of  re¬ 
fracting  and  diffracting  objects  from  data  obtained  in  a 
manner  analogous  to  ordinary  computed  tomography. 
Emphasis  has  been  on  cases  where  time  of  flight  or 
optical  path  length  data  are  recorded  for  radiation 
passing  through  a  continuous  refracting  object.  Be¬ 
cause  this  problem  is  nonlinear  it  has  been  approached 
by  iterative  schemes  or  by  linearized  inverse  scattering 
theory.  Only  a  few  empirical  studies  of  these  ap¬ 
proaches  have  been  carried  out  to  date.  Although 
some  of  these  studies  have  been  successful,  there  is 
much  to  be  learned  about  their  convergence  behavior, 
reliability,  etc.  The  subject  is  challenging,  interesting, 
and  of  potential  importance  in  applications  ranging 
from  plasma  diagnostics  to  medical  ultrasonic  imag¬ 
ing. 

Readers  should  note  that  I  have  not  discussed  close¬ 
ly  related  topics  such  as  electromagnetic  and  seismic 
geophysical  exploration-’'23  or  addressed  the  related 
problem  of  refractive  errors  in  ultrasonic  tomography 
based  on  attenuation  measurements. 

I  would  like  to  thank  Richard  Gordon  for  suggesting 
this  review  and  for  helpful  comments.  This  work  was 
sponsored  in  part  bv  the  U.S.  Army  Research  Office. 
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This  paper  was  presented  at  the  symposium  on  In¬ 
dustrial  Applications  of  Computed  Tomography  and 
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Tomographic  reconstruction  of  a  function  from  the  values 
of  its  line  integrals  is  a  well-known  diagnostic  technique  for 
use  in  various  fields.1  Severn!  reconstruction  algorithms 
have  been  developed  with  the  assumption  that  the  probing 
radiation  propagates  along  straight  lines.  However,  if  re¬ 
fractive  index  gradients  normal  to  the  direction  of  wave 
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Fig.  1.  Schematic  representation  of  the  formation  of  an  interfero- 
gram  with  appreciable  ray  bending.  Hay  DFFF.  first  exposure, 
taken  with  the  index  of  refraction  nn  uniform  everywhere.  Rav 
ARCP':  second  exposure,  taken  with  an  object  of  unknown  index  of 
refraction  n(x,v)  varying  inside  a  circular  zone  of  radius  rn.  Imaging 
system  focuses  on  plane  o  -  o. 


propagation  cause  significant  ray  bending,  the  use  of  straight 
line  inversion  algorithms  may  lead  to  unsatisfactory  recon¬ 
structions.  Ways  to  correct  for  refraction  effects  have  been 
the  subject  of  active  research,2  9  and,  at  present,  this  pheno- 
mennm  may  he  taken  into  account  through  the  iterative  use 
of  lengthy  digital  rav  tracing  algorithms.2"9  Alternatively,  a 
perturbation  technique  is  available  which  does  not  require 
the  use  of  computational  ray  tracing.  This  approach  was 
originally  developed  by  Norton  and  Linzer9  for  use  in  ultra¬ 
sonic  imaging.  In  this  Letter  we  present  a  modification  of 
that  technique,  applicable  to  holographic  and  Mach- 
Zehnder  interferometry  of  transparent  objects. 

Figure  1  is  a  schematic  representation  of  the  formation  of 
an  interferogram  by  two  initially  parallel  plane  waves.  It  is 
assumed  that  outside  the  circular  region,  and  throughout 
during  one  of  the  exposures,  the  index  of  refraction  has  the 
uniform  value  tin.  The  undeviated  ray  DF.FP',  representa¬ 
tive  of  that  exposure,  interferes  with  ray  ARCP’.  which  goes 
through  the  object.  The  lens  in  Fig.  1  is  focused  on  the  object 
plane  a  -  o.  Point  P  in  that  plane  is  the  apparent  origin  of 
both  rays.  Point  f’  is  defined  on  the  refracted  ray  so  that  its 
x  coordinate  is  equal  to  that  of  the  unrefracted  ray  at  the 
point  where  it  leaves  the  test  section.  J^oint  Fjs_defined  on 
the  unrelracted  rav  so  that  distances  PC  and  PF  are  equal. 
The  optical  pathlength  difference  (OPD)  between  the  two 
rays  is  given  by 


AT  =  nda  -  n„(xf  -  xn )  , 


(11 


w  here  n  is  the  unknown  index  of  refraction  of  the  object  and 
where,  if  C  is  a  generic  point  in  Fig.  1,  x<;  denotes  its  x 
coordinate.  Suppose  there  was  no  refraction.  Then  the 
OPI)  would  be 


AT  =  (n  -  ri„ldx  . 


(2) 


where  F  i«  the  point  at  which  both  rav s  enter  the  test  section. 
Following  Norton  and  Linzer.'1  we  now  determine  a  relation 
between  AT  and  AT. 

Assume  that  the  index  of  refraction  deviates  only  slightly 
from  the  ambient: 


.'  A  V  • 


n(x,y)  «  n0  +  <h(x,y)  , 


where  <  is  a  small  parameter.  Assume  also  that  the  curved 
trajectory  can  be  expressed  as  a  small  perturbation  of  a 
straight  line: 


'  »|1  +  »/(x) 


(3) 


where  yn  is  the  y  coordinate  of  point  P.  If  the  function  htx.y) 
is  expanded  about  the  line  y  =  yp  and  if  we  use  Eq.  (3),  an 
approximation  to  the  value  of  the  index  of  refraction  along 
the  ray  trajectory  is  obtained.  Substitution  of  this  approxi¬ 
mation  into  Eq.  (1)  and  use  of  Eq.  (2)  yield,  after  some 
rearrangement,9 


AT  “  AT  —  c  , 


where 


XC  1 

r 2yp(/hJ>  +  7yntfpf'2)dl  ~  no(*F  ~  xc> 

*r  ^ 


and  where  h denotes  the  partial  derivative  of  h(x,y )  with 
respect  to  y  evaluated  at  y  =  yp. 

Norton  and  Linzer9  obtained  an  expression  for  f(x)  by 
using  Fermat's  principle  to  minimize  the  integral  rids  and 
using  standard  techniques  from  the  calculus  of  variations. 
The  result  is 


f(x)  - 


-Mr 

Wp  [he 


(x  -  p)htyp)dii  +  a,x  +  a,  . 


The  integration  constant  a!  is  easily  evaluated  by  noting  that 
f'{x)  =  0  forx  <  xp,  so  that  a)  =  0.  To  determine  the  value  of 
the  constant  a2  the  following  geometric  condition  is  used: 


dy 

dx 


yc  -  yp 


,m’c  xc  -  Xp 


After  some  algebra,  the  perturbation  of  the  trajectory  is 
found  to  be 


/<*>' 


— —  ("f  (x  -  p)hf,(p)dp  -  I  (xP-  p)h$(p)dp\,  (4) 

ndyp  [h,  I*,  J 

which  differs  from  the  expression  developed  by  Norton  and 
Linzer  for  ultrasonic  imaging  because  different  boundary 
conditions  apply.  Equation  (4)  can  now  be  used  to  evaluate 
the  term  c.  Noting  from  Fig.  1  that 


xr-xr  =  (x, 


■  -  x„)  — - - 1 

[costT) 


that  for  small  T 


1 


cos(</>) 


- 1 


-  tan  (</>)  ' 


IO’f/'(*c>|2 


and  defining 

H(x)  =  n0[o(/'(x)|2  ’ 


p  12 

I  ihp{p)dp 


we  obtain 


1 

C  =  2 


(xr  -  xr)H{xr)  ■ 


H(x)dx 


\ 


(5) 


The  computation  of  this  correction  term  is  based  on  a  first 
estimate  of  n(x.y)  obtained  by  straight  line  inversion  of  AT. 
After  r  is  computed,  the  modified  OPDs  are  determined  and 
straight  line  inverted  to  yield  an  improved  estimate  for  the 
index  of  refraction. 

Consider,  as  an  example,  a  refractive-index  distribution  of 
the  form 
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Fig.  2.  Root -mean  square  error  of  the  uneorrected  (solid  line)  and 
corrected  (dashed  line)  straight  line  reconstructions  vs  object  plane 
selected  by  the  imaging  system.  The  location  of  the  object  plane 
relative  to  the  center  of  the  field  is  constant  for  each  experiment. 
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Fig.  3.  Cross  section  of  the  test  field  through  line  x  =0.  Solid  line: 
actual  profile.  Dashed  line:  straight  line  reconstructed  profile. 
Focusing  plane  at  ip  =  0.0. 


I 


Fig  4.  Cross  section  of  the  test  lield  through  line  t  =  0  Solid  line: 
actual  profile.  Dashed  line  perturbation  corrected  profile.  Same 
focusing  plane  as  in  Fig.  3. 


n 

"n 


1 .0  -  0.!).r>  exp  — 


1 1-  +  tv  -  ivirj 
0.09 


0.05  exp  — 


[x‘  +  (\  +  0.5F| 
0.04 


in  which  the  coordinate  ip.its  are  arhitrnrv.  The  OPDs  were 
simulated  hv  numerically  inlegraling  the  rav  equation  to 
find  the  actual  rav  traiectories  through  this  field.  We  traced 


thirty-two  rays  per  view  for  each  of  thirty-two  viewing  direc¬ 
tions.  It  was  found  that  with  r0  =  1.0  the  maximum  bending 
angle  is  8.5°.  To  use  Eq.  (5),  xp  needs  to  be  specified.  We 
varied  this  distance  from  -r0  to  +r0  and  kept  it  constant  for 
each  computational  experiment.  For  all  reconstructions  we 
used  the  well-known  filtered-backprojection  algorithm10  on 
a  40  X  40  grid.  The  root  mean  square  error  of  the  reconstruc¬ 
tion  is  defined  as 


_ 100  /T  y 

maxjn^,  -  nj  y  m  Z 


I",,  ~  ".>) 


where  n,,  is  the  exact  value  of  the  field  in  pixel  i,  nir  is  the 
reconstructed  value  in  that  same  pixel,  and  m  is  the  total 
number  of  pixels.  The  error  resulting  from  reconstructions 
ignoring  refraction  is  the  solid  line  in  Fig.  2.  No  results  are 
shown  for  xp  <  — 0.2ro  as  ray  crossing  occurred  beyond  that 
position  of  the  object  plane.  After  the  perturbation  correc¬ 
tion  for  refraction  is  applied,  the  error,  as  shown  by  the 
dashed  line  in  Fig.  2,  decreases  significantly. 

Another  way  to  visualize  the  performance  of  this  method  is 
to  compare  a  cross  section  of  the  reconstructed  refractive- 
index  profiles.  This  has  been  done  in  Figs.  3  and  4,  for  which 
we  have  chosen  xp  =  0.  In  those  figures  the  actual  profile  is 
the  solid  line,  while  the  dashed  line  is  the  reconstructed 
profile.  Figure  3  shows  the  result  of  the  reconstruction 
ignoring  refraction.  It  can  be  seen  that  the  corrected  recon¬ 
struction  shown  in  Fig.  4  is  much  closer  to  the  actual  profile. 

If  the  index  of  refraction  varies  along  one  space  coordinate 
only,  as  in  a  boundary  layer,  the  above  equations  reduce  to  a 
much  simpler  form.  Details  and  an  example  for  this  case  will 
appear  in  a  forthcoming  paper. 

This  work  was  presented  at  the  OSA  Annual  Meeting,  Wash¬ 
ington,  DC,  Oct.  1985.  This  research  is  sponsored  by  the 
Army  Research  Office. 
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Refraction  correction  in  holographic  interferometry  and 
tomography  of  transparent  objects 


Ignacio  H.  Lira  and  Charles  M.  Vest 


In  this  paper  we  review  and  extend  the  state  of  the  an  in  the  algorithms  that  have  been  developed  to 
tomographically  reconstruct  1-D  and  2-D  refractive-index  fields  in  the  presence  of  significant  refraction.  A 
perturbation  approach  and  two  iterative  procedures  were  tested  and  compared  in  numerical  simulation  of 
holographic  interferometry  experiments.  Due  to  the  nonlinearity  of  the  problem,  it  is  very  difficult  to  draw 
general  conclusions  with  respect  to  the  behavior  of  the  iterative  algorithms,  which  is  divergent  in  the  examples 
presented  here.  In  contrast,  the  perturbation  technique,  which  is  the  easiest  one  to  implement  and  the  fastest 
to  run,  is  shown  to  be  very  powerful  in  reducing  refraction  errors. 


I.  Introduction 

Inverting  the  data  acquired  through  holographic  in¬ 
terferometry  to  reconstruct  refractive- index  fields  is  a 
well-established  diagnostic  technique.1  By  and  large, 
the  underlying  approximation  in  the  development  of 
the  conventional  reconstruction  algorithms  is  that  the 
probing  radiation  travels  along  a  straight  line,  chang¬ 
ing  in  phase  but  not  in  direction,  as  it  crosses  the 
refractive-index  inhomogeneities.  In  reality,  when 
light  encounters  a  refractive-index  gradient  normal  to 
its  direction  of  propagation,  it  is  refracted.  There  are 
many  situations  in  which,  for  practical  purposes,  ray 
bending  may  be  neglected.  However,  ignoring  refrac¬ 
tion  often  produces  appreciable  reconstruction  errors. 
This  has  been  shown  to  be  the  case  not  only  in  interfer¬ 
ometry2  but  also  in  related  areas  such  as  geophysical3 
and  ultrasonic4  5  imaging. 

Analytical  solutions  to  the  highly  nonlinear  problem 
of  tomographic  reconstruction  of  strongly  refracting 
objects  have  not  been  found,  so  one  must  resort  to 
computational  techniques.  In  a  tutorial  paper,  Vest2 
discussed  several  aspects  of  the  problem  and  outlined 
some  of  the  approaches  which  have  been  taken  to  find 
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an  approximate  solution.  In  this  paper  a  more  com¬ 
prehensive  and  detailed  review  is  presented.  In  Sec.  II 
we  address  the  1-D  case,  which  has  been  considered 
extensively  by  previous  investigators.  Section  III  dis¬ 
cusses  the  axisymmetric  problem,  which  has  also  re¬ 
ceived  considerable  attention.  In  Sec.  IV  the  more 
general  2-D  asymmetric  problem  is  considered.  Sec¬ 
tion  V  contains  the  results  of  numerical  experiments 
designed  to  investigate  and  compare  the  performance 
of  the  refraction  correction  schemes.  Section  VI  sum¬ 
marizes  the  results  obtained. 

II.  One-Dimensional  Refractive-Index  Fields 

We  show  in  Fig.  1  the  geometry  of  the  data  collection 
arrangement  for  the  analysis  of  a  1-D  boundary  layer 
using  holographic  interferometry  with  plane  wave  illu¬ 
mination.  Above  a  solid  object  of  length  L  there  is  a 
region  with  a  refractive-index  variation  of  the  form 
n(y)  for  y  <  5.  It  is  assumed  that  outside  the  boundary 
layer,  and  throughout  the  entire  region  during  one  of 
the  exposures,  the  index  of  refraction  has  the  uniform 
value  n<).  The  undeviated  ray  DPCP,  representative 
of  that  exposure,  interferes  with  a  refracted  ray  ABP. 
The  lens  in  Fig.  1  is  focused  on  the  object  plane  o-o. 
Point  P  in  that  plane  is  the  apparent  origin  of  both 
rays.  Point  B  is  the  point  at  which  the  refracted  ray 
leaves  the  test  section.  Point  C  is  defined  on  the 
unrefracted  rav  such  that  the  distance  FT  is  equal  to 
the  distance  Pb.  Due  to  the  presence  of  the  lens,  the 
optical  path  lengths  from  BtoP  and  from  C  to  P  are 
equal.  The  optical  path  length  difference  <OPD)  be¬ 
tween  the  two  rays  is,  then. 
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and  where  x c  is  the  x  coordinate  of  point  C.  The 
subscript  in  the  term  <J>p  emphasizes  the  fact  that,  to 
the  observer,  a  fringe  located  at  P  appears  to  be  pro¬ 
duced  by  the  interference  of  two  overlapping  straight 
rays  through  point  P.  The  ray  path  AB  is  obtained 
from  the  ray  equation6 

which,  in  2-D  Cartesian  coordinates,  reduces  to 

av"(S'£y)11+yi'  <3) 

where  y  =  y(x)  is  the  path  of  the  ray  and  the  super¬ 
scripts  indicate  derivative  with  respect  to  x.  Using  the 
initial  condition  y  =0  for  x  =  0  and  assuming  n  =  n(y), 
the  integration  of  this  equation  gives 


where  nA  is  the  index  of  refraction  at  y  =  y.-,.  Since  the 
ray  path  element  is  ds  -  yl  +  y'1,  the  optical  path 
length  of  the  refracted  ray  through  the  test  section 
can  be  expressed  as  either 

(l+/2)dx,  (5) 


1  fL  -w 
i  —  n-ax. 

"a  Jo 


One  of  the  earliest  papers  to  assess  and  correct  for 
the  influence  of  refraction  was  that  of  Wachtel.7  That 
work  was  subsequently  extended  and  modified  by  the 
research  of  Howes  and  Buchele.6-11  Other  papers  in 
which  this  situation  has  been  considered  are  Refs.  12- 
17. 

Howes  and  Buchele11  based  their  analysis  on  a  Tay¬ 
lor  series  expansion  of  the  index  of  refraction  around 
its  value  at  the  point  where  the  rays  enter  the  test 
section,  i.e., 

n  m  nA +  y  f>ity  -  (6) 

The  general  analysis11  is  quite  involved  and  ne>— 
be  repeated  here.  The  equations  are  very  simple 
ever,  if  only  two  terms  in  expansion  (6)  are  kept,  t  is, 

if 

n  -  nA  +  n\(y  -yA),  (7) 

where  n\  is  the  derivative  of  the  index  of  refra  ion 
evaluated  at  y  =  v  A.  The  ray  equation  has,  in  this  case, 
the  following  exact  solution: 


y  -y^  +  i^cosh^^-  i  ■  (8) 

where  lA  3  nA/n'A.  If  it  can  be  assumed  that  L/lA  «  1, 
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Fig.  1.  Interferometric  analysis  of  a  1-D  boundary  layer-type  re¬ 
fractive-index  field. 


expanding  the  hyperbolic  cosine  and  keeping  only  the 
first  two  terms  in  the  expansion,  Eq.  (8)  becomes 

1  X2 

ymy*  +  lTA-  (9) 

To  this  first  level  of  approximation,  then,  the  ray 
traces  are  parabolas.  Differentiating  Eq.  (9)  and  sub¬ 
stituting  the  result  into  Eq.  (5)  give  the  following  ap¬ 
proximation  for  the  optical  path  length  of  the  refracted 
ray: 


■H©]- 


On  the  other  hand,  from  the  geometry  of  Fig.  1  it  is 
readily  seen  that  the  optical  path  length  of  the  unre¬ 
fracted  ray  is 

*o  ”  "o*c  ”  "o(*P  +  <£•  -  */>W 1  +  tan*0|.  (11) 

where  <t>  is  the  exit  angle  of  the  ray.  Using  Eq.  (9),  we 
can  approximate  the  tangent  to  this  angle  by 

tan*  at  j-  ■  (12) 


Putting  this  expression  into  Eq.  (11),  expanding  the 
square  root,  and  substituting  the  result,  together  with 
Eq.  (10),  into  Eq.  (1),  we  obtain 


where  we  have  defined  K  =  1  -  xp/L. 

The  index  of  refraction  at  y  -  yA  may  now  be  ob¬ 
tained  from  Eq.  (13)  provided  n'A  is  known.  Howes 
and  Buchele  suggest  that  n'A  can  be  found  by  differen¬ 
tiating  Eq .  ( 1 3 )  and  assuming  that  the  derivative  of  the 
second  term  on  the  right-hand  side  is  small  compared 
nth  that  of  the  first  term.  Doing  this  we  find 

l  da* I  .... 

n'.  at-— — |  (14) 

L  dy  |>*.v 

so  that  the  refractive-index  gradient  is  proportional  to 
the  measurable  fringe-shift  gradient.  Alternatively, 
the  effect  of  the  second  term  can  be  accounted  for  by 
using  Eq.  (14)  as  a  first  approximation  to  n'A,  substi¬ 
tuting  this  value  into  Eq.  (13),  and  differentiating 
again  to  obtain  a  refined  estimate.  This  process  can  be 
repeated  if  necessary. 

The  problem  with  this  method  is  that  the  actual 
entrance  coordinates  are  unknown  to  the  observer.  If 
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Fig.  2.  Schematic  representation  of  the  formation  of  an  mterfero- 
gram  in  a  2-D  axisvmmetric  refractive-index  field 


one  uses  coordinate  y p,  the  resulting  distribution  be¬ 
comes  distorted.  Recognizing  this  problem  Howes 
and  Buchele8  showed  that  within  the  linear  approxi¬ 
mation  the  difference  yp  —  y  disappears  when  K  = 
1/2,  i.e.,  when  the  imaging  system  is  focused  at  the 
center  of  the  test  section.  This  can  easily  be  seen  to  be 
the  case.  From  Fig.  1  we  have 

\'p  -  Vp 

tarn#  * - -  (15) 

Taking  xp  =  L  and  using  the  parabolic  approximation. 
Eqs.  (9)  and  1 12),  we  find 


which  verifies  the  assertion. 

In  summary,  if  the  linear  approximation,  Eq.  (7),  is 
valid  for  the  range  y.s  <  yp,  Howes  and  Buchele  recom¬ 
mend  using  Eq.  (13)  together  with  Eq.  (14)  (with  itera¬ 
tion,  if  necessary)  to  find  n(y).  The  data  should  be 
obtained  with  the  imaging  system  focused  at  the  center 
of  the  test  section. 

Higher-order  refraction  correction  is  possible  if 
more  terms  in  Eq.  (6)  are  included  and  higher-order 
derivatives  of  the  fringe  order  number  are  calculated. 
Because  the  fringes  are  obtained  experimentally,  the 
accuracy  of  these  derivatives  becomes  increasingly  du¬ 
bious. 

A  different  approach  to  account  for  refraction  was 
taken  by  Svensson.18  In  his  analysis  he  concluded 
that,  with  the  proper  choice  of  the  object  plane,  the 
fringe  data  could  be  inverted  satisfactorily  as  if  no 
refraction  had  taken  place.  Vest1  simplified  Svens- 
son’s  analysis,  but  his  derivation  contains  an  error  in 
an  intermediate  step.  The  correct  equations  are  as 
follows.  Let  A4>  be  the  actual  OPD  and  let  A$  be  the 
OPD  that  would  be  observed  if  no  refraction  had  taken 
place.  Clearly, 

A<frp  *  I  rip  —  n0iL,  (17) 

but,  from  Eq.  (7), 

npasn,  +  n  ,( \>  —  _v4i.  1 18) 

Combining  Eqs.  1 17)  and  ( 13)  and  using  Eqs.  ( 18)  and 
(16),  we  find 


Thus,  if  it  can  be  assumed  that  nJnA  ~  1  throughout 
the  boundary  layer,  the  difference  between  A4>  and  A$ 
is  minimized  when  K  =  1/3,  i.e.,  when  the  object  plane 
is  chosen  to  be  at  a  distance  of  1/3  of  the  section  length 
from  the  window  closest  to  the  imaging  system.  In 
other  words,  at  that  position  of  the  object  plane  the 
fringe  count  may  be  straight  line  inverted,  as  if  no 
refraction  had  occurred.  Note  that,  contrary  to 
Howes's19  assertion,  minimizing  refraction  errors  in 
this  way  does  not  result  in  a  distorted  distribution  of 
the  index  of  refraction.  Empirical  confirmation  of 
this  theory  has  been  given  in  Ref.  15. 

III.  Axially  Symmetric  Refractive-Index  Fields 

There  are  many  instances  of  practical  importance  in 
which  the  index  of  refraction  may  be  considered  essen¬ 
tially  constant  along  one  direction,  while  having  axial 
symmetry  in  a  plane  perpendicular  to  that  direction. 
This  occurs,  for  example,  when  forming  interfero- 
grams  of  density  in  flames,  thermal  jets,  and  plumes,  or 
of  plasma  electron  density  around  exploding  wires,  arc 
discharges,  or  laser  beams. 

Figure  2  shows  a  circular  test  section  of  radius  r0 
being  probed  by  two  plane  waves.  One  of  them,  with 
representative  ray  DPC.  corresponds  to  the  exposure 
with  uniform  index  of  refraction  n,i  everywhere.  The 
other  wave,  with  representative  ray  AB,  corresponds  to 
the  exposure  taken  with  the  object  present,  for  which 
the  index  of  refraction  is  of  the  form  n(r).  We  assume 
that  during  both  exposures  the  index  of  refraction  in 
the  medium  outside  of  the  test  section  has  the  value  n0. 
The  imaging  system — not  shown  in  the  figure — is  fo¬ 
cused  on  an  object  plane  o-o  located  a  distance  x  P  from 
the  center  of  the  field.  Point  P  is  the  conjugate  of  the 
point  in  the  image  plane  at  which  both  rays  interfere. 

Kahi  and  Mylin20  wrote  the  OPD  of  these  rays  in  the 
following  form: 

«  Qp  +  2 n0zA  +  tan*  +  n„xF (sec®  -  1),  (20) 

where  4>p  =  n(r)ds  and  where  the  refraction  angle  <t> 
must  be  taken  as  positive  if  the  ray  bends  upward,  and 
negative  otherwise.  The  first  two  terms  on  the  right- 
hand  side  of  this  equation  represent  the  OPD  that 
would  be  obtained  if  no  refraction  had  taken  place,  the 
third  term  is  a  contribution  due  to  refraction,  and  the 
fourth  term  is  an  additional  focusing  effect.  Interest¬ 
ingly,  this  equation  shows  that  all  refraction  contribu¬ 
tions  are  known  if  <j>  is  known.  The  OPD  given  by  Eq. 
(20)  is  a  function  of  the  yP  coordinate  given  by 

yp  »  yA  seep  +  xp  tar i«.  (21) 

According  to  Kahl  and  Mylin,  if  a  conventional  Abel 
inversion  algorithm  is  used  to  reconstruct  the  refrac¬ 
tive-index  field,  it  is  better  to  focus  on  the  center  of  the 
test  section,  so  that  the  third  contribution  in  Eq.  (20)  is 
eliminated  and  the  difference  between  the  refracted 
and  unrefracted  OPDs  is  kept  small.  Vest21  conduct¬ 
ed  an  extensive  computer  simulation  of  interferomet¬ 
ric  experiments  using  several  analytic  functions  n(r). 
All  the  functions  studied  had  zero  slope  at  r  =  0  and 
smoothly  approached  n0  at  r  =  r0.  He  found  that  Abel 
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inversion  of  the  actual  interferometric  data  A4>/'  pro¬ 
duces  remarkably  accurate  reconstructions  when  ip  = 
0.  The  above  argument  explains  in  part  this  result. 
Vest's  findings  were  later  confirmed  in  Refs.  22-24. 

In  Ref.  21  it  is  shown  that 

■t>p  -  2  p - f - ;  dr.  (22) 

"p  r y  ij-  -  p- 

where  n  =  rn(r),  p  =  and  rP  is  obtained  by  solving 
v(rp)  —  p.  In  the  same  reference  it  is  also  shown  that 
the  an  ie  between  OA  and  OB  subtended  by  the  re¬ 
fracted  -ay  is  given  by 

h-2  p— f—  dr.  (23) 

‘'pry/t r  -  p: 

Using  the  Abel  transform.  Vest21  inverted  Eqs.  (22) 
and  (23).  He  found  that  the  index  of  refraction  could 
be  obtained  from  either 

or 

If"0  ,  i’p\  da  j  1 

rinl  =  r  exp  -j  cnsh  (  j— <ip  •  '24) 

provided  r(  17)  is  a  single-valued  function  of  radius. 

It  is  easy  to  show  that,  if  the  exponent  in  Eq.  (24)  is 
integrated  by  parts,  the  following  simpler  alternative 
result  is  obtained: 

r(,)-r0eip[-i  P-^p=^pl  •  (25) 

L  VP'  -  J 

The  index  of  refraction  can  then  be  obtained  from  this 
equation  if  the  central  angle  0  is  known.  This  result 
was  derived  by  Maruyama  et  al .25  and  independently 
by  Zimin  and  Frik.26  The  former  authors  recommend 
the  following  procedure  to  find  0.  First,  obtain  the 
refraction  angle  <t>  from  the  derivative  of  the  fringe 
pattern: 

c/A'I’d 

sind  *  -  • 

(J  V  p 

Next,  calculate  the  incident  angle  in  by  means  of  the 
following  geometric  relation: 


vertical  copper  cylinder.  Comparison  with  thermo¬ 
couple  measurements  gave  excellent  agreement. 

A  similar  reconstruction  procedure  was  derived  by 
Hunter  and  Schreiber.28  They  showed  that  the  index 
of  refraction  can  be  obtained  directly  from  the  fringe- 
shift  data  by  using  the  following  equation: 

1  C*  . !p\  1  dA<t>"  .  1 

nlr|)“n.jexp - I  cosh  -) - - — up  ■  '26) 

[  «■  J„  \ljp  dp 

Note  that  the  data  are  taken  as  a  function  of  yP  but  the 
derivative  is  with  respect  to  v,*.  However,  Hunter  and 
Schreiber  tacitly  assumed  that  the  apparent  and  actu¬ 
al  entrance  coordinates  are  equal.  F rom  the  geometry 
of  Fig.  2  it  may  be  seen  that  this  situation  will  occur  for 
a  certain  location  of  the  object  plane,  but  this  location 
is  not  the  same  for  all  ray  pairs.  The  applicability  of 
Eq.  (26)  is  therefore  restricted  to  the  case  in  which 
refraction  is  very  small. 

IV.  Asymmetric  2-D  Fields 

We  turn  now  to  the  more  general  case  of  2-D  asym¬ 
metric  refractive-index  fields.  This  problem  has  re¬ 
ceived  little  attention  in  interferometry,  the  only 
works  known  to  us  being  the  iterative  algorithm  oi  Cha 
and  Vest-2  and  the  perturbation  approach  of  Lira  and 
Vest.10  Additionally,  another  iterative  approach  has 
been  proposed  in  the  areas  of  ultrasound  tomogra¬ 
phy5-31  and  in  geophysics3-32  but,  to  our  knowledge,  it 
has  not  been  used  in  interferometry.  The  techniques 
are  as  follows. 

A.  Straight  Line  Inversion  with  Modified  Data  (SUM) 

This  procedure  was  first  introduced  in  ultrasound 
tomography  by  Bates  and  McKinnon33  34  and  in  inter¬ 
ferometry  by  Cha  and  Vest.29  Let  us  first  define  a4>  as 
the  actual  OPD  including  refraction,  obtained  by  nu¬ 
merical  or  physical  experiment,  and  A4>  as  the  OPD 
that  would  be  obtained  in  a  hypothetical  experiment 
involving  the  same  object,  but  where  no  refraction 
occurs.  The  subscript  P  has  been  dropped  for  concise¬ 
ness.  Cha  and  Vest’s  procedure  consists  of  an  iterative 
scheme  in  which  a  series  of  estimates  of  A#  are  ob¬ 
tained,  hopefully  converging  to  its  exact  unknown  val¬ 
ue.  The  algorithm  is 


Finally,  obtain  the  central  angle  from  0  =  ir  -  <f>  -  2i0. 
Maruyama  et  al.  successfully  tested  this  procedure  in  a 
simulated  experiment  in  which  n(r)  was  given  analyti¬ 
cally25  and  with  an  acrylic  homogeneous  cylinder  im¬ 
mersed  in  a  fluid  with  a  slight  mismatch  in  the  index  of 
refraction.27 

Zimin  and  Frik.26  manipulating  Eq.  (25),  put  it  in  a 
form  more  convenient  for  calculations.  Their  result  is 


They  used  this  equation  to  measure  the  temperature 
profile  in  the  thermal  boundary  layer  around  a  heated 


A<t*+I  =  A<f*  +  (A4  -  A4*)  with  k  =  1,2 .  (27) 

where  A$*  is  the  kth  estimate  of  A<$,  n*  is  the  refrac¬ 
tive-index  distribution  resulting  from  the  straight  line 
inversion  of  A$*.  and  Al>*  is  the  OPD  obtained  compu¬ 
tationally  through  digital  ray  tracing35  routines  over 
the  field  n*. 

It  is  usually  convenient  to  start  by  letting  Al*1  *  A<$. 
i.e.,  refraction  is  initially  ignored  and  the  data  are 
straight  line  inverted  to  yield  n 1 .  Cha  and  V est29  were 
very  successful  in  applying  this  technique  to  1-D  re¬ 
fractive-index  fields.  However,  in  two  dimensions  the 
results  were  less  conclusive.  It  was  found  from  nu¬ 
merical  simulations  that  this  algorithm  will  tend  to 
decrease  the  reconstruction  error  in  the  first  few  itera¬ 
tions,  but  eventually  it  diverges. 
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B.  Perturbation  Approach 

A  similar,  but  computationally  much  simpler  tech¬ 
nique,  has  been  presented  by  Norton  and  Linzer0  for 
ultrasonic  applications.  In  their  approach  it  is  as¬ 
sumed  that  the  index  of  refraction  differs  only  slightly 
from  its  value  at  the  surroundings  and  that  the  ray 
trajectory  deviates  very  little  from  the  straight  line. 
With  these  assumptions  an  improved  estimate  of  Ad*  is 
obtained  and  subsequently  straight  line  inverted, 
without  the  need  for  the  computational  expense  of  ray 
tracing  and  iteration.  We30  have  modified  Norton  and 
Linzer's  analysis  due  to  differences  between  the  data 
acquisition  systems  used  in  ultrasonic  imaging  and  in 
holographic  interferometry.  Our  result  is 

A$  at  A‘£  —  c. 


where 


xf,>mxr) 


H{x)dx 


and  where 


HU)  5  - 


n^>(  )du 


In  this  expression  n>,  is  the  derivative  of  n(x,y)  with 
respect  to  v  evaluated  at  y  =  yp.  Points  C,  E,  and  Pare 
defined  in  Fig.  3.  The  computation  of  H(x)  must  be 
based  on  a  first  estimate  of  n(x,y).  As  in  the  SLIM 
algorithm,  this  estimate  may  be  obtained  by  straight 
line  inversion  of  the  actual  OPD  data  A<i. 

If  the  index  of  refraction  varies  along,  say,  the  y 
coordinate  only,  the  above  equations  reduce  to  a  much 
simpler  form.  Using  the  points  shown  in  Fig.  1,  it  is 
easy  to  show07  that 


(28) 


Note  that  in  this  case  c  becomes  zero  if  the  object  plane 
is  chosen  such  that  xp  =  2xp/3.  This  result  supports 
the  discussion  following  Eq.  (19). 


C.  Curved  Ray  Algebraic  Inversion  (CRAI) 

A  second  iterative  technique  to  obtain  the  index  of 
refraction,  also  in  the  context  of  ultrasound,  was  first 
suggested  by  Johnson  et  al.5  and  later  by  Schomberg.01 
The  basis  of  this  procedure  is  the  utilization  of  the 
current  estimate  of  the  index  of  refraction,  nk,  to  ob¬ 
tain  a  system  of  algebraic  equations,  where  the  un¬ 
knowns  are  the  corrections  for  this  estimate  at  each 
pixel.  Since  this  technique  has  not  been  used  in  inter¬ 
ferometry,  we  proceed  to  describe  the  algorithm  in 
more  detail. 

The  optical  path  length  for  a  digital  ray  inside  the 
test  section  is  approximated  by 

K,-l 

■t"  -fir  ,  I  -1-  h  V  n‘(r +  hnk(TlK  ),  (29) 

i-2 

where  points  r,,  with  l  =  1.2, ...JC,  are  obtained  by 
digital  ray  tracing i:’  over  the  field  nK.  In  this  equation 
subscript  i  identifies  a  particular  ray  and  h  is  the 
distance  between  the  points.  The  length  h ,  depends 


Fig.  3.  Schematic  representation  of  the  formation  of  an  interfero- 
gram  in  a  2  D  asymmetric  refractive-index  field. 


on  the  location  of  point  Tk,  with  respect  to  the  circle 
boundary.  Since  in  general  the  points  r,i  will  not 
coincide  with  the  pixel  centers,  interpolation  is  needed 
to  find  the  values  of  the  index  of  refraction  at  those 
points.  In  general,  we  can  write 

m 

n*(r,,)  =  V  cjjrt*  l  =  1.2 _ JC..  (30) 

;TT 

where  n*  is  the  value  of  the  index  of  refraction  at  pixel  j 
in  the  sth  reconstruction  and  m  is  the  number  of 
pixels.  The  coefficients  c*;  depend  on  the  interpola¬ 
tion  scheme.  F or  example,  with  bilinear  interpolation 
only  the  coefficients  corresponding  to  the  four  nearest 
pixels  are  nonzero.  Substitution  of  Eq.  (30)  into  Eq. 
(29)  gives 

m 

*?  “  £  a‘>n''  (31) 
/-> 

where 

K,-‘ 

a‘  “  T'0  +  H  1  +  (3Z) 

.-2 

The  rav  whose  optical  path  is  given  by  Eq.  (29)  inter¬ 
feres  with  a  straight  ray  for  which  an  appropriate  dis¬ 
tance  (f  must  be  found,  such  that 

+0.  =  not,*.  (33) 

We  now  form  the  difference  between  the  optical  path 
lengths  given  by  Eqs.  (31)  and  (33)  for  several  rays  in 
each  viewing  direction  and  arrange  the  result  as  a 
vector  which  we  will  call  Ad>*.  Thus, 

A**  -  A*n*  -  n„l*.  (34) 

where  A*  is  a  matrix  whose  components  are  given  by 
Eq.  (32),  n'1  is  a  vector  formed  by  the  nk  pixel  values, 
and  1*  is  a  vector  whose  components  are  the  distances 

£ 

Equation  (34)  suggests  that  the  actual  OPDs  can  be 
written  as 


A<f  •  An  -  n0l, 


(35) 
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where  n  is  the  vector  of  exact  pixel  values.  Matrix  A 
and  vector  1  are  unknown  quantities. 

Let  us  define  a  correction  for  the  index  of  refraction 
as  follows:  (in*  s  n  -  n*.  We  now  substitute  this 
definition  into  Eq.  t35>  and  assume  that  A  A*  and 
that  1  as  1*.  If  the  result  is  combined  with  Eq.  (34)  we 
obtain 


A*  in* 


A4>  -  A4>* 


1 36) 


This  system  of  equations  must  be  solved  for  the  correc¬ 
tion  6n*.  after  which  the  next  estimate  of  the  index  of 
refraction  is  computed  according  to  n**1  =  n*  +  5n*. 

We  conclude  this  section  by  noting  that,  since  Eq. 
(36)  represents  a  system  which  will  usually  be  large, 
sparse,  and  either  under  or  over  determined  (depend¬ 
ing  on  the  number  of  rays  traced  and  on  the  grid  size), 
the  solution  should  be  obtained  with  a  series  expansion 
method.38 

V.  Numerical  Experiments 

In  this  section  the  behavior  of  the  refraction  correc¬ 
tion  algorithms  is  compared  and  evaluated  by  means  of 
computational  experiments  involving  two  simple  trial 
functions.  Additional  results  are  presented  in  Ref.  37. 

A.  Two-Dimensional  Example 
Consider  the  following  refractive-index  field: 

x2 -Hy  -  0.1  )2] 


(37) 


mx^y)  =  n.,  -  0.01no  (exp  - 


+  exp  - 


”  {  r  0.09 

[x2  +  (y  +  0.5)2]  I 
0.04  }  ' 


The  OPD  data  for  this  double  Gaussian  distribution 
was  generated  using  the  ray  tracing  procedure  in  Ref. 
37,  with  fifty  rays  for  each  of  twenty  viewing  directions, 
spaced  9°  apart.  If  the  interferogram  is  formed  using 
an  imaging  system  focused  at  the  center  of  the  test 
section,  there  will  be  about  thirty-six  fringes  across  the 
object  using  a  He-Ne  laser  and  with  the  millimeter  as 
the  unit  of  length.  This  fringe  number  corresponds  to 
a  view  in  the  direction  of  the  y  axis.  It  was  found  that 
the  maximum  bending  angle  at  the  exit  of  the  refrac¬ 
tive-index  field  is  1.73°. 

In  the  following  study  we  assumed  that  the  imaging 
system  was  kept  fixed  while  the  object  is  rotated  about 
the  center,  so  that  the  same  value  of  xp  results  for  all 
viewing  directions. 

Figure  4  shows  the  straight  line  reconstructed  re¬ 
fractive-index  profile  for  xp  =  2.0.  In  this  figure  the 
ordinate  axis  represents  the  difference  1000  X  (n  -  rc0) 
in  a  cross  section  through  the  plane  x  =  0,  in  which  the 
solid  line  is  the  exact  profile  and  the  dashed  line  is  the 
reconstruction.  The  reconstruction  error  was  mea¬ 
sured  in  terms  of  both  the  rms  and  maximum  error 
percentages,  defined  as 


el  rms)  = 


100  /  1  v  .  " 

- - - \/  —  N  tn.,-n.rr, 
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Fig.  4.  Straight  line  reconstructed  profile  (dashed  line)  of  the  dou¬ 
ble  Gaussian  distribution  compared  to  the  actual  profile  (solid  line). 
The  focusing  plane  is  at  xp  =  2.0.  In  obtaining  this  figure,  twenty 
viewing  angles  were  considered,  with  fifty  rays  per  view.  Resolution 
is  40  x  40  pixels.  The  root-mean  square  error  of  this  reconstruction 
is  1.76%  and  the  maximum  error  is  8.74%. 
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Fig.  5.  Perturbation-corrected  reconstruction  idashed  line)  of  the 
double  Gaussian  distribution  compared  to  the  actual  profile  (solid 
iine).  Root-mean  square  error  is  0.35%.  maximum  error  is  2.97%. 


where  n,,.  is  the  exact  value  of  the  field  at  the  center  ol 
pixel  i,  n,r  is  the  reconstructed  value  at  that  same  pixel, 
and  m  is  the  number  of  pixels.  F or  this  reconstruction 
the  rms  error  is  1.76%  and  the  maximum  error  is  8.74%. 

I.  Results  with  the  Perturbation  Approach 

After  the  perturbation  correction  is  applied,  the  re¬ 
constructed  profile  looks  much  closer  to  the  actual 
profile,  as  can  be  observed  in  Fig.  5.  In  the  corrected 
profile  the  rms  and  maximum  errors  are  0.35%  and 
2.97%,  respectively,  which  means  a  factor  of  5  in  error 
reduction.  In  obtaining  these  figures,  both  straight 
line  inversions  were  done  with  the  filtered  backprojec- 
tion  algorithm39  with  a  Shepp- Logan40  filter  on  a  grid 
of  40  X  40  pixels. 

The  quality  of  the  uncorrected  reconstruction  varies 
quite  significantly  as  the  object  plane  changes.  Figure 
6  shows  the  rms  and  maximum  errors  as  a  function  of 
the  object  plane.  In  this  figure  the  solid  lines  give  the 
errors  obtained  for  the  first  (uncorrected)  reconstruc- 
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Fig.  6.  Reconstruction  error  as  a  function  of  the  object  plane  select¬ 
ed  for  the  straight  line  (solid  lines)  and  perturbation -corrected 
idaahed  lines i  reconstructions  of  the  double  Gaussian  function. 
Plain  line  give  rms  errors,  while  lines  with  asterisks  give  maximum 
errors. 
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Fig.  8.  Reconstruction  error  as  a  function  of  inversion  number  for 
the  reconstruction  of  the  double  Gaussian  function  using  the  CRAJ. 
The  solid  line  gives  the  rms  error,  while  the  dashed  line  shows  the 
maiimum  error.  F or  comparison,  the  perturbation-corrected  errors 
ior  the  second  inversion  are  shown  as  white  and  black  symbols  for 
the  maiimum  and  rms.  respectively.  The  focusing  plane  is  ip  «  2.0. 
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Fig.  7.  Reconstruction  error  as  a  function  of  inversion  number  for 
the  reconstruction  of  the  double  Gaussian  function  using  the  SLIM 
iterative  approach  The  solid  line  gives  the  rms  er-or,  while  the 
lashed  line  shows  the  maximum  error.  For  comparison,  the  pertur- 
lation-correcied  errors  tor  the  second  inversion  are  shown  as  white 
ind  black  symbols  for  the  maximum  and  rms.  respectively  The 
focuaing  plane  is  xo  =  2.0. 


tion,  while  the  dashed  lines  are  the  errors  for  the  inver¬ 
sion  obtained  after  the  correction  is  applied.  The 
plain  lines  give  rms  errors,  while  the  lines  with  aster- 
iska  are  maximum  errors.  Results  are  shown  only  for 
-2.5  <  xp<  7.5,  as  ray  crossing  occurred  beyond  those 
positions  of  the  object  plane.  From  this  Figure  it  may 
be  seen  that  the  error  caused  by  refraction  is  almost 
negligible  when  the  object  plane  is  slightly  in  front  of 
the  center  of  the  Field,  and  that  the  uncorrected  enors 
are  very  sensitive  to  the  location  of  the  object  plane. 
In  fact,  these  errors  are  very  nearly  directly  propor¬ 
tional  to  the  object  plane  distance  from  the  plane 
where  refraction  error  is  minimized.  Thus,  to  mini¬ 
mize  the  refraction  error  the  imaging  system  has  to  be 
focused  very  carefully  close  to  the  center  of  the  distri¬ 
bution.  The  perturbation  corrected  errors  are  not  as 


sensitive  to  the  object  plane  location.  As  Fig.  6  shows, 
they  are  almost  constant  over  the  range  —2<xp<2. 

2.  Results  with  the  SLIM 

The  SLIM  procedure  yields  the  results  in  Fig.  7. 
Since  this  procedure  is  iterative,  it  is  convenient  to 
show  its  behavior  by  plotting  the  rms  and  maximum 
error  percentages  as  a  function  of  the  inversion  num¬ 
ber  (where  inversion  1  is  the  first,  uncorrected  recon¬ 
struction)  for  a  given  object  plane.  Space  limitations 
preclude  the  presentation  of  equivalent  plots  for  other 
objects  planes.  We  believe,  however,  that  one  object 
plane  should  be  sufficient  to  acquire  a  good  idea  of  the 
main  trends  in  the  behavior  of  this  algorithm.  In  Fig.  7 
the  focusing  plane  has  been  selected  to  be  xp  =  2.0. 
For  direct  comparison  with  the  perturbation  ap¬ 
proach,  all  results  have  been  obtained  using  twenty 
views,  fifty  rays  per  view  and  40  X  40  pixels.  In  this 
example  the  SLIM  reduces  the  reconstruction  error  in 
the  first  few  iterations,  but  after  inversion  5  the  proce¬ 
dure  diverges. 

3.  Results  with  the  CRA1 

In  the  following  results  we  have  retained  the  number 
of  views,  the  number  of  rays  per  view,  the  grid  size,  and 
the  focusing  plane  so  that  direct  comparison  between 
the  three  approaches  can  be  done.  The  CRAI  algo¬ 
rithm  has  some  additional  degrees  of  freedom,  namely, 
the  type  of  algebraic  reconstruction  procedure,  the 
number  of  iterations  in  that  procedure,  and  the  relax¬ 
ation  factor  to  be  used  (which  may  vary  from  one 
iteration  to  the  next).  In  this  example  the  reconstruc¬ 
tions  were  obtained  using  the  SART  algorithm  of  An¬ 
dersen  and  Kak41  with  two  iterations  ( i.e.,  the  pixels 
were  updated  at  the  end  of  the  second  SART  iteration) 
and  a  constant  relaxation  factor  of  0.5. 

As  shown  in  Fig.  8,  the  CRAI  algorithm  yields  a 
mildly  unstable  solution  when  reconstructing  the  dou- 
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Fij?.  9.  Temperature  distribution  in  the  boundary  layer  surround¬ 
ing  a  heated  vertical  Hat  piate  submerged  in  water.  Plate  tempera¬ 
ture  is  29°C  and  ambient  water  temperature  is  20°C.  Plate  width  is 
10  cm  and  the  results  shown  are  tor  a  plane  i  cm  above  the  leading 
edge.  In  this  figure  the  object  plane  is  at  the  point  where  the  actual 
and  apparent  plate  boundaries  coincide.  The  solid  curve  shows  the 
temperature  as  given  in  Ref.  42,  using  which  the  interferometric  data 
was  generated.  The  dasned  curve  shows  the  temperature  obtained 
with  the  conventional  straight  line  inversion  of  the  data.  The  curve 
with  asterisks  gives  the  perturbation-corrected  temperature. 


ble  Gaussian  function,  diverging  after  approximately 
the  tenth  iteration. 

B.  One-Dimensional  Example 

In  this  subsection  the  performance  of  the  refraction 
correction  schemes  for  the  case  of  1-D  ray  bending  is 
investigated.  The  chosen  example  simulates  the  re¬ 
fractive-index  field  resulting  from  temperature  varia¬ 
tions  in  a  natural  convective  boundary  layer  about  a 
vertical  flat  plate  at  29°C  submerged  in  water  at  20°C. 
The  analytic  form  of  this  field  is  given  in  detail  in  Ref. 
37  and  need  not  be  repeated  here.  Our  ray  tracing 
model37  predicts  a  maximum  bending  angle  of  2.13°  for 
the  innermost  ray.  In  this  example  we  have  chosen  the 
object  plane  such  that  the  apparent  and  actual  loca¬ 
tions  of  the  wall  coincide. 

1.  Results  with  the  Perturbation  Approach 

The  result  of  subtracting  the  correction  term  as  giv¬ 
en  by  Eq.  r281  from  the  OPD  data  is  shown  in  Fig.  9.  In 
this  figure,  the  solid  curve  represents  the  exact  tem¬ 
perature  profile,  the  dashed  curve  represents  the  tem¬ 
perature  obtained  using  the  conventional  straight  line 
inversion  of  the  OPD  data,  and  the  curve  with  asterisks 
gives  the  perturbation-corrected  inversion.  There  is  a 
very  noticeable  improvement  as  a  result  of  applying 
the  correction  term. 

2.  Results  with  the  SLIM 

In  implementing  the  iterative  algorithms  to  reduce 
refraction  errors  in  the  1-D  case,  the  discrete  form  of 
the  ray  tracing  routines  require  the  current  estimate  of 
the  index  of  refraction  in  the  form  of  a  pixel  array  in 
which  each  entry  is  the  value  of  the  index  of  refraction 
at  a  given  y  coordinate.  Here,  and  in  the  next  subsec¬ 
tion,  50  pixels  have  been  used.  Figure  10  shows  the 
rms  error  and  the  absolute  value  of  the  maximum  error 
for  each  straight  line  inversion  of  the  corrected  data. 
The  perturbation -corrected  errors  are  also  shown  in 
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Fig.  10.  Reconstruction  error  as  a  function  of  inversion  number  for 
the  reconstruction  of  the  1-D  temperature  profile  using  the  SLIM 
approach.  The  solid  line  gives  the  rms  error,  while  the  dashed  line 
shows  the  maximum  error.  For  comparison,  the  perturbation  cor- 
rected  errors  for  the  second  inversion  are  shown  as  white  and  black 
symbols  for  the  maximum  and  rms.  respectively.  Percentages  are 
based  on  50  pixel  values.  All  other  parameters  are  as  in  Fig.  9. 
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Fig.  11.  Reconstruction  error  as  a  function  of  inversion  number  for 
the  reconstruction  of  the  I  D  temperature  protile  using  the  ORAL 
The  solid  line  gives  the  rms  error,  while  the  dashed  line  shows  the 
maximum  error  F or  comparison .  the  perturbation-corrected  errors 
for  the  second  inversion  are  shown  as  white  and  black  symbols  for 
the  maximum  and  rms.  respectively.  Percentages  are  based  on  50 
pixel  values.  All  other  parameters  are  as  in  Fig.  9. 


that  figure,  as  was  done  in  sec.  V.  A.  It  is  seen  in  Fig.  10 
that  the  SLIM  procedure  produces  a  substantial  error 
decrease  up  to  inversion  number  eight,  but  afterward 
the  algorithm  diverges. 

3-  Results  with  the  CRAI 
As  shown  in  Fig.  11,  the  convergence  of  the  CRAI  for 
this  example  is  very  fast  compared  to  the  SLIM  proce¬ 
dure.  Thus,  at  the  fourth  inversion  the  errors  are 
0.22%  rms  and  0.93%  maximum,  and  the  resulting  pro¬ 
file  practically  coincides  with  the  actual  one.  A  small 
undulation  in  the  reconstructed  profile,  however,  leads 
to  computational  ray  crossing  when  we  attempt  to 
continue  the  iterations  after  the  fourth  inversion. 
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VI.  Conclusions 

Due  to  the  nonlinearity  of  the  problem,  it  is  very 
difficult  to  draw  general  conclusions.  F rom  the  exam¬ 
ples  presented  above,  and  from  other  examples  dis¬ 
cussed  in  Ref.  37,  it  is  seen  that  the  quality  of  the 
reconstructions  and  the  convergence  and  divergence 
characteristics  of  each  algorithm  are  strongly  object- 
dependent.  The  most  interesting  result  of  this  re¬ 
search  is  the  fact  that  the  perturbation  technique, 
which  is  the  easiest  one  to  implement,  is  quite  powerful 
in  reducing  refraction  errors.  This  approach  has  the 
important  advantage  of  not  requiring  the  use  of  digital 
ray  tracing,  thereby  significantly  reducing  computa¬ 
tional  costs  compared  to  the  iterative  procedures. 
Moreover,  it  is  very  appealing  since  it  makes  use  of 
straight  line  inversion,  a  technique  supported  by  an 
extensive  body  of  knowledge. 

The  SLIM  algorithm  is  essentially  a  less  sophisticat¬ 
ed  attempt  to  find  the  difference  between  actual  and 
unrefracted  OPDs.  According  to  Eq.  (27),  the  condi¬ 
tion  for  the  convergence  of  this  algorithm  is  that  the 
difference  A<f>  —  must  grow  progressively  smaller. 
The  examples  presented  here  have  shown,  however, 
that  in  some  cases  that  assumption  is  not  correct.  To 
explain  this  behavior,  consider  the  following  argu¬ 
ment.  Suppose  that  eventually  the  exact  solution  is 
found,  i.e.,  that  for  some  k  we  have  nh  =  n.  Because  of 
discretization  and  round-off  errors,  the  calculated 
OPDs  A4>*  will  differ  from  the  exact  OPDs  A4>  by  a 
finite  amount.  The  next  estimate  of  the  OPDs  with¬ 
out  refraction,  A4*+1,  will  then  be  different  from  the 
ones  that  produced  the  correct  reconstruction.  As  a 
consequence  the  resulting  index  of  refraction,  n*+i, 
will  be  different  from  n.  Since  this  difference  is  not  a 
product  of  refraction,  it  is  not  likely  to  be  reduced  after 
further  iterations. 

The  CRAI  technique  for  reconstruction  along 
curved  rays  is  based  on  a  totally  different  approach  to 
the  problem.  It  does  not  use  straight  line  inversion 
techniques,  and  therefore  it  is  more  difficult  to  imple¬ 
ment  and  slower  than  the  SLIM.  Nevertheless,  we 
have  showed  that  the  CRAI  is  a  very  effective  tech¬ 
nique  in  correcting  for  refraction  errors.  According  to 
Eq.  (36),  the  requirements  for  the  convergence  of  this 
algorithm  are  the  same  as  that  for  the  SLIM.  Thus, 
because  numerical  errors  compound  after  the  mini¬ 
mum  error  solution  is  found,  the  CRAI  will  also  tend  to 
diverge  after  an  acceptable  solution  has  been  found. 
The  divergence  of  the  CRAI  is,  however,  usually  less 
pronounced  than  that  obtained  with  the  SLIM. 

In  conclusion,  the  evidence  presented  here  suggests 
that,  if  refraction  becomes  non-negligible.  the  pertur¬ 
bation  correction  should  be  applied.  For  comparison, 
it  may  be  convenient  to  use  one  of  the  iterative  ap¬ 
proaches.  in  which  case  the  CRAI  algorithm  is  recom¬ 
mended.  We  feel,  however,  that  further  research  is 
needed  to  establish  the  conditions  in  which  the  itera¬ 
tive  algorithms  will  converge  or  diverge. 
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Abstract.  Digital  holographic  interferometry  is  a  hybrid  optical- 
digital  technique  for  determining  the  phase  of  an  interferogram. 
This  technique  improves  the  accuracy  of  interferometric  measure¬ 
ment  of  fluid  properties  and  enhances  the  utility  of  interferom¬ 
etric  flow  visualization.  Displays  of  the  interferometric  phase 
produce  excellent  images  of  weakly  refracting  two-dimensional 
flows  and  can  be  used  to  produce  integral  projection  images  of 
three  dimensional  flows  which  differ  from  and  -ompiement 
schlieren  and  shadowgraph  images.  The  technique  is  explained 
herein  and  examples  of  its  use  in  both  continuous  wave  and 
pulsed  interferometry  are  presented. 


1  Introduction 

Schlieren,  shadowgraph  and  interferometric  flow  visual¬ 
ization  techniques  use  the  integrated  effect  of  a  fluid 
optical  property  on  a  beam  of  light  passing  through  an 
object  to  form  an  image  of  a  flow  pattern.  These  tech¬ 
niques  work  well  for  a  variety  of  flows  but  have  limited 
utility  for  three-dimensional  flows  and  weakly  refracting 
flows.  Schlieren  and  shadowgraph  images  are  formed  by 
ray  bending  which  is  approximately  proportional  to  re¬ 
fractive  index  gradients  and  second  derivatives  of  the 
refractive  index,  respectively.  The  latter  two  techniques 
work  well  for  two-dimensional  flows  characterized  by 
large  refractive  index  gradients,  such  as  shock  patterns. 
They  also  have  been  used  to  observe  complex  structures 
like  turbulent  jets.  Such  images  show  qualitative  flow 
features  such  as  the  gross  outline  of  the  motion,  the  fined 
grained  structure  of  turbulence  (Crow  &  Champagne 
1971),  and  the  presence  of  large,  two  dimensional  struc¬ 
tures  [e  g.  the  vortex-like  structures  of  Brown  &  Roshko 
(1974)].  but  conceal  other  structural  aspects  of  the  flow 
such  as  the  presence  of  unmixed,  entrained  fluid  in¬ 
clusions  within  the  flow. 

Interferometric  methods  are  based  on  the  phase  delay 
of  a  plane  wave  passing  through  the  object  and  produce 
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an  integral  projection  of  the  object's  refractive  index  field 
However,  this  phase  delay  is  encoded  in  a  fringe  pattern 
which  can  be  difficult  to  interpret.  The  fm  -es  are  useful 
for  visualizing  objects  such  as  plumes  and  boundary 
layers,  but  become  complicated  for  three-dimensional 
objects.  None  of  these  methods  works  well  for  weakly- 
refracting  flows  because  the  gradients  involved  produce 
indistinct  schlieren  and  shadowgraph  images  and  yield 
broad,  ambiguous  interferometric  fringes. 

Digital  interferometry  is  a  technique  by  which  the 
interferometric  phase  delay  is  determined  very  accurately 
at  a  large  number  of  points  in  the  image.  The  phase  may 
be  displayed  as  a  gray  scale,  yielding  an  image  that  has 
none  of  the  ambiguities  associated  with  Schlieren  and 
shadowgraph  images.  Furthermore,  the  technique  may  be 
used  to  image  weakly-refracting  flows  since  direct  deter¬ 
mination  of  the  phase  allows  one  to  enhance  the  contrast 
of  flow  details  that  can  not  be  resolved  by  conventional 
interferometry.  Finally,  the  technique  may  be  used  to 
make  distributed  measurements  of  flow  quantities. 

2  Digital  interferometry 

Digital  interferometry'  is  a  recently  developed  hybrid 
optical-digital  metrology  technique  combining  two  expo¬ 
sure  holographic  interferometry  with  digital  image 
acquisition  and  computer  processing  to  determine  the 
interferometric  phase  directly  from  a  set  of  image  irradi- 
ance  measurements  (Dandliker  &  Thalmann  1985:  Flan- 
haran  1985).  This  technique  is  similar  to  heterodyne 
holographic  interferometry  in  that  both  manipulate  the 
interferogram  s  phase  in  a  known  manner  to  determine  its 
magnitude.  The  image  intensity  of  a  holographic  inter- 
ferogram  is  given  by  (Vest  1979): 

/  (  v.  i  l  =  /n(  v.  r)  ;  1  +  m  (  v,  i  )  cos[<£  (  a.  j  )  «*][  (  I  1 

where  /,, ( v.  r >  is  the  background  intensity,  m(.v.  i  )  is  the 
fringe  contrast  and  is  the  interferometric  phase. 

The  term  w  is  the  uniform  phase  bias  term  which  m  the 
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case  of  digital  or  heterodyne  interferometry  can  be  mani¬ 
pulated.  The  image  irradiance  is  recorded  for  several 
different  values  of  this  bias  phase  using  a  digitizing 
camera  and  is  stored  in  computer  memory. 

The  unknown  phase  <t>  can  be  calculated  from  the 
ulues  of  ihe  recorded  image  irradiance  distribution  and 
ihe  known  values  of  <p  using  any  one  of  a  number  of 
algebraic  relationships  reviewed  recently  by  Hariharan 
1 1985)  in  this  study,  we  recorded  n  separate  irradiance 
patterns  with  the  reference  phase  if  evenly  distributed 
from  zero  to  2 The  phase  <t>  may  then  be  determined, 
modulo  2  r.  hv : 

n 

^  /.  (  v.  v )  sin  if, 

i  v.  i  i  =  tan'  - -  (2) 

T  /  (  v.  v)  cos  if, 

- 1 

w  nere 

/ix.w  =  /  i.v.  y )  '  1  ■+•  rn  i.x.  v)  cos  [0  (v.  i)  +  <j>\\  (3) 

and  o  =  — — .  n  =  number  of  exposures.  The  inverse  tan- 
n 

cent  function  maps  the  fringe  pattern  into  a  linear,  dis¬ 
continuous  function.  This  enables  one  to  eliminate  the 
phase  Mgn  ambiguity  normally  associated  with  cosinuso¬ 
idal  fringes.  The  usual  fringe  counting  procedure  is 
•eplaced  by  a  simple,  computational  sorting  operation. 
A  negative  discontinuity  indicates  an  increase  in  fringe 
number  and  a  positive  discontinuity  a  decrease.  Because 
the  phase  is  evaluated  independently  at  each  point  in  the 
mage,  its  determination  is  unaffected  by  spatial  variation 
n  the  background  irradiance  or  fringe  contrast.  The 
accuracy  of  phase  determination  may  be  of  the  order  of 
!  50-  I  100  of  a  fringe,  compared  with  1/100- 1/1.000  of 
a  fringe  for  heterodyne  interferometry  and  1/5-1/10 
'-.nge  for  conventional  interferometry,  depending  in  ail 
cases  on  the  nature  of  the  object  being  studied  (Dandliker 
&  Thalmann  1985).  The  combined  use  of  computational 
fringe  counting  and  high  resolution  image  storage  devices 
enables  one  to  resolve  complex  fringe  patterns.  Digital 
nterferometry  is  a  convenient,  accurate  high  resolution 
interferometric  technique. 

This  method  can  be  used  for  both  real-time  (Hariharan 
et  al  1982)  and  double  exposure  holographic  mterferorn- 
ctrv  Because  real-time  images  show  temporal  phase 
variation  due  to  even  small  perturbations  which  are 
nherent  in  many  fluid  mechanics  experiments,  double 
exposure  methods  are  preferred  for  most  steady  flow 
Judies  and  all  unsteady  flow  studies.  The  double  expo- 
-ure  technique  requires  two  reference  waves,  one  for  each 
exposure  (Fig.  1 1.  The  first  exposure  is  made  without  the 
hiect  i flow  field)  and  the  second  exposure  is  made  with 
■he  object  present.  The  film  is  developed  and  the  image  is 
•eoonstructed  by  illuminating  it  with  both  reference  waves 
omultaneousiy  Primary,  conjugate  and  cross  reconstruc¬ 
ts  are  then  present  (Dandliker  el  al.  I97fi|.  The  con¬ 


jugate  and  cross  reconstructions  must  be  properly  com¬ 
pensated  for  so  that  their  presence  does  not  significantly 
reduce  the  accuracy  of  the  technique  (Dandliker  et  al. 
1982).  The  two  primary  reconstructions  overlap  to  form 
the  desired  interferogram  with  the  phase  bias  term  given 
by: 

if  =  k  (r,  —  r;)  -  k  (/•;  -  r’:)  (4) 

2  ?i 

where  k  =  —:  /.  =  recording  wavelength;  rt.r.  are  the 
/. 

reference  source  distances  during  recording,  and  r\.  r:  are 
the  reference  source  distances  during  reconstruction.  The 
phase  bias  term  if  may  be  shifted  by  changing  the  path 
length  of  either  reference  wave  by  a  small  amount.  Jr 
Therefore  the  change  in  the  phase  bias  term  is  J  if  =  k  J  r 

This  technique  can  be  applied  to  flow  visualization  or 
measurement  with  either  plane-wave  or  diffused  object 
illumination.  The  plane  wave  setup  (Fig.  2)  is  preferable 
when  laser  power  .s  limited.  When  a  plane  object  wave  is 
used,  the  two  reference  waves  must  have  a  wide  angular 
separation  in  order  to  eliminate  troublesome  overlap  of 
the  cross  reconstructions. 

Diffused  illumination  (Fig.  3)  is  generally  preferable, 
however,  because  it  allows  the  two  reference  waves  to 
have  a  small  angular  separauon.  thereby  reducing  the 
errors  due  to  misalignment  of  the  hologram  with  the  two 
reference  waves.  .Although  the  cross-reconstructions  over¬ 
lap.  they  may  be  almost  completely  decorrelated  if  their 
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digital  interferogram 


images  are  displaced  in  the  image  plane  by  an  amount 
greater  than  the  speckle  size.  This  overlap  reduces  the 
fringe  contrast,  but  does  not  significantly  reduce  the 
accuracy  iDandliker  et  al.  1982).  The  small  angular 
separation  of  the  two  reference  waves  causes  a  regular, 
visible  fringe  pattern  to  be  present  on  the  hologram 
surface,  which  causes  a  periodic  error.  This  effect  is 
minimized  by  placing  a  low  f-number  imaging  lens  very 
close  to  the  hologram  so  that  the  regular  fringe  partem  is 
outside  the  field  of  focus. 

Breuckm3nn  and  Thieme  (1985)  used  this  arrangement 
to  correct  for  errors  due  to  wavelength  shift  when  per¬ 
forming  digital  interferometry  when  the  recording  and 
reconstruction  lasers  had  different  wavelengths  (see 
Fig.  4a  and  b).  In  this  case,  the  bias  phase  is  given  by: 


i>  -  (A|  r:-  k:r'z)  *  k:.)r: 
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2  n  2  n 

where  kt  =  — ,  ki  =  — ,  /.t  =  recording  wavelength.  = 

f.2 

reconstruction  wavelength,  or 

.Up  =  k]  (r,  -  r:)  —  Aw(r|  —  rt)  +  A;  Jrt . 

For  reference  beams  separated  by  a  small  angle  <5 
whose  bisector  forms  an  angle  0  w  ith  the  hologram  plane: 

r,  -  /■-  =  r  cos  6  sin  <3 . 


cos  0  sin  <3,  cos  sin  3- 


_  ,  ,  ,  cost?  sin <3 

The  terms  of  the  form - ; - can  be  shown  to  be  the 

frequency  of  the  fringe  pattern  formed  b\  the  two  refer¬ 
ence  waves  in  the  hologram  plane.  Therefore,  by  setting 
such  that  the  fnnge  spacing  of  the  reconstruction  refer¬ 
ence  waves  is  equal  to  a  pattern  that  wouid  be  formed  by 
the  two  recording  reference  waves,  the  chromatic  errors  in 
the  phase  are  eliminated.  This  allows  the  interferometric 
phase  to  be  evaluated  in  the  usual  wav 


3  Experimental  setup 

A  digital  interferometer  with  a  plane  object  wave  is  shown 
in  Fig.  2.  The  two  reference  waves  are  separated  by  about 
20°.  The  hologram  was  held  and  developed  in  a  real-time 
liquid  gate  to  reduce  alignment  errors.  A  second  set  of 
experiments  was  made  with  diffused  illumination  (Fig.  4) 
using  a  pulsed  ruby  laser  for  recording  and  a  He-Ne  laser 
for  reconstruction  in  a  manner  described  by  Breuckmann 
and  Thieme  (1985).  The  second  recording  reference  wave 
was  derived  in  this  case  by  tilting  a  mirror  in  the  original 
reference  beam  by  about  0.5°.  The  hologram  is  recon¬ 
structed  with  a  He-Ne  laser  using  a  Michelson  interferom¬ 
eter  to  create  the  desired  fringe  pattern.  In  both  cases  the 
phase  shifting  is  done  by  translating  a  mirror  in  one  of  the 
reference  beams  a  fraction  of  a  wavelength.  The  mirror  is 
mounted  on  a  piezoelectric  cell,  to  which  a  control  voltage 
is  applied.  In  both  cases,  the  two  reference  waves  are 
brought  together  to  form  a  fringe  pattern  which  serves  as 
a  measure  of  their  mutual  phase  difference.  A  pair  of 
photodiodes  are  placed  in  the  fringe  pattern  and  the  error 
signal  between  them  is  used  with  a  proportional-integral 
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controller  to  generate  the  control  voltage  for  the  piezo¬ 
electric  cell  to  maintain  a  stable  fringe  pattern  (and 
therefore  phase  shift)  in  the  presence  of  electrical,  me¬ 
chanical  and  thermal  perturbations.  The  plane  wave  inter- 
ferographic  images  were  made  with  a  128  by  128  pixel 
array  CID  camera  while  the  difTuse-illumination  images 
were  made  with  a  video  camera  in  conjunction  with  a  256 
by  384  pixel  wdeo  frame  store.  Each  pixel  has  255  gray 
levels.  Image  computations  were  done  on  an  LSI-1 1  based 
computer,  which  is  interfaced  with  the  camera. 

4  Results 

The  plane-waxe  interferometer  was  used  to  visualize  the 
flow  around  a  cylindrical  heated  wire.  Figure  5a  show's 
the  digitized  interferogram  of  a  strong,  laminar  plume 
while  Fig.  5b  shows  the  gray  scale  display  of  the  phase, 
modulo  2  n.  The  conversion  of  the  fringe  pattern  to  a 
linear  discontinuous  function  is  evident.  Figure  5c  shows 
the  digitized  interferogram  of  a  slightly  heated  wire  in  a 
cross-flow.  This  interferogram  consists  of  a  few  narrowly 
spaced  fnnges  near  the  wire  and  a  single  broad,  indistinct 
fringe  in  the  wake.  The  phase  display  of  this  inter¬ 
ferogram  (Fig.  5 d)  has  a  greatly  enhanced  contrast  that 
clearly  shows  thermal  variation  in  the  wake,  which  is  not 
apparent  in  the  interferogram. 

The  diffuse-illumination  interferometer  was  used  to 
visualize  a  nominally  axisymmetnc.  turbulent,  helium  jet 
injected  in  still  air.  The  interferogram  shown  in  Fig.  6  a  is 
complicated  and  difficult  to  interpret  especially  near  the 
jet  center  where  the  fringes  become  broad  and  indistinct. 
The  phase  display  (Fig.  6b)  is  an  improvement,  producing 
a  contour  map  of  the  phase.  Using  a  sorting  procedure 
based  on  the  nature  of  discontinuities  in  the  phase 


display,  we  can  compute  the  total  phase  at  each  point  in 
the  image.  Displaying  this  total  phase  as  a  gray  scale 
(normalized  to  255  gray  levels)  yields  the  image  in  Fig.  6c. 
This  image  represents  a  true  integral  projection  of  the 
density  of  the  jet.  Close  examination  of  this  image  reseals 
two  interesting  features.  The  first  is  the  axial  distribution 
of  areas  of  high  concentration  of  helium.  The  second  ib 
the  presence  of  large  inclusions  of  unmixed  ambient  fluid 
near  the  center  line  of  the  jet.  Although  the  implications 
of  these  structures  will  not  be  discussed  here,  it  is 
interesting  to  compare  the  absolute  phase  plot  to  the 
schlieren  images  of  a  similar  jet  made  by  Crow  and 
Champagne  (1971)  and  the  laser-induced  fluorescence 
images  of  Dimotakis  et  al.  (1983).  The  schlieren  pictures 
show  an  apparent  fine-grained  "surface”  of  the  jet  with 
little  indication  of  the  internal  structure.  The  integral 
phase  plot  is  similar  to  the  centerline  laser  fluorescence 
pictures,  in  that  both  indicate  areas  of  high  jet  fluid  con¬ 
centration  and  large  inclusions  of  unmixed  ambient  fluid 
This  is  somewhat  surprising  because  the  integration  would 
tend  to  average  out  the  effects  of  locally  strong  variation 
in  the  concentration  field.  Further  study  of  such  images 
could  be  used  to  examine  the  symmetry  properties  of  the 
jets  and  other  mixing  phenomena. 

5  Applications 

These  examples  point  out  the  potential  utility  of  digital 
interferometry  for  flow  visualization  and  measurement.  It 
greatly  enhances  the  image  contrast  of  interferograms  of 
weakly  refracting  flows,  removes  the  sign  ambiguity  as¬ 
sociated  with  conventional  interferograms.  and  produces 
images  that  display  rather  subtly  flow  features.  These 
features  could  allow  unique  applications  of  interferomeir\ 
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to  flow  visualization,  such  as  using  small  temperature 
variations  as  tracers  of  vorticity  in  non-buoyant.  two- 
dimensional  flows  and  generation  of  stereoscopic  images 
of  three-dimensional  scalar  fields.  Furthermore,  the  tech¬ 
nique's  accuracy,  high  spatial  resolution,  and  facility  of 
operation  indicate  potential  for  making  spatially  distrib¬ 
uted  scalar  measurements. 

Two-dimensional  scalar  fields  may  be  evaluated  direct- 
l>  from  the  phase,  the  experimental  geometry  and  the 
constitutive  relation  of  the  refractive  index  to  the  scalar 
quantity  of  interest  (eg.  the  Gladstone-Dale  relation). 
Cross-sections  of  three-dimensional  fields  may  be  evalu¬ 
ated  using  optical  tomography,  a  technique  for  recovering 
the  cross-sectional  distributions  of  a  field  from  its  integral 
projections  (Sweeney  &  Vest  1973.  Vest  1979).  A  large 
number  of  contiguous  cross-sections  could  be  obtained 
from  a  set  of  projection  images,  allowing  three  dimen¬ 
sional  maps  of  the  field  to  be  calculated.  This  would  allow 
investigation  of  the  structure  of  the  scalar  distribution. 
We  have  a  system  for  doing  this  underdevelopment. 

The  technique  also  has  the  potential  for  determining 
the  concentration  in  multicomponent  mixtures  through 
the  use  of  simultaneous  mterferograms  with  several  dif¬ 
ferent  wavelengths.  Since  each  component  has  a  different 
specific  refractivity  at  each  wavelength,  the  phase  mea¬ 
surement  in  each  wavelength  is  a  linear  combination  of 
the  concentration  of  each  constituent,  implying  a  unique 
relation  between  the  phase  measurements  and  the  consti¬ 
tuent  concentrations.  The  technique's  sensitivity  and  its 
ability  to  register  efficiently  the  various  mterferograms 
could  make  such  a  system  practical. 
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